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ABSTRACT 

We use Laskar's frequency mapping technique to study the dynamics of triaxial galaxies 
with central density cusps and nuclear black holes. For ensembles of ~ 10 4 orbits, 
we numerically compute the three fundamental frequencies of the motion, allowing us 
to map out the Arnold web. We also compute diffusion rates of stochastic orbits in 
frequency space. The objects of greatest importance in structuring phase space are 
found to be the 3-dimensional resonant tori, regions where the fundamental frequencies 
satisfy a relation of the form = Ilo\ + mioi + nuj^ with integer (I, m, n). When stable, 
resonant tori generate phase space regions in which the motion is regular; these regions 
are not necessarily associated with a stable periodic orbit as in systems with only 
two degrees of freedom. Boxlike orbits are generically stochastic but some tube orbits 
are stochastic as well. The spectrum of diffusion rates for boxlike orbits at a given 
energy is well approximated as a power law over at least six decades. Models with high 
central concentrations - steep central cusps or massive black holes - exhibit the most 
stochasticity. Even a modest black hole, with a mass ~ 0.3% the mass of the galaxy, is 
as effective as the steepest central density cusp at inducing stochastic diffusion. There 
is a transition to global stochasticity in boxlike phase space when the mass of a central 
black hole exceeds ~ 2% of the galaxy mass. We estimate the dependence of orbital 
evolution rates on galaxy structural parameters. We predict a greater average degree 
of dynamical evolution in faint elliptical galaxies due to their high central densities 
and short crossing times. The evolution time is estimated to be shorter than a galaxy 
lifetime for absolute magnitudes fainter than about —19 or —20, consistent with the 
observed change in many elliptical galaxy properties at this luminosity. 



1 INTRODUCTION 



Motion in a smooth gravitational field becomes quite sim- 
ple if the number of isolating integrals equals or exceeds the 
number of degrees of freedom, and much work in galactic dy- 
namics has focussed on finding integrable or near-integrable 
models 



The large core of the Perfect Ellipsoid is a poor match to 
real elliptical galaxies, all of which ex hibit power-law cusps 
in the stellar density at s mall radii ( | Fcr rarosc_ et al. 1994 ; 
Merritt fc Fridman 199E|; |Gebhardt ct al. 1996|). There is 



increasingly strong evidence that many elliptical galaxies 
and bulges also contain massive central objects, possibly the 



fui galactic polenLials. Kuzmiu (1950, 1973) showed 
that th ere is a unique, ellipgoidally-stratificd mass model for 

which the corresponding potential has three global integrals 
of the motion, quadratic in the velocities. Kuzmin's model 
- explored in detail by de Zeeuw (1985) who christened it 
the "Perfect Ellipsoid" - has a large, constant-density core 
in which the orbital motion is that of a 3-D harmonic os- 
cillator. Every orbit in the core of the Perfect Ellipsoid fills 
a rectangular parallelepiped, or box. At higher energies in 
the Perfect Ellipsoid, box orbits persist and three new orbit 
families appear: the tubes, orbits that preserve the direction 
of their circulation around either the long or short axis of the 
figure. Tube orbits respect an integral of the motion anal- 
ogous to the angular momentum, and hence - unlike box 
orbits - avoid the center. All in tegrable triax ial potentials 
have a similar orbital structure ( Hunter 1995 ). 



black holes that are though t to have powered quasars (Ko- 



rmendy & Richstone 1995). While the masses of these dark 



central components are often very uncertain, typical esti- 
mates are 1(T 3 < M h /M g < 1(T 2 , where M h is the black 
hole mass inferred from the orbital motions of surrounding 
stars and gas and M g is the stellar mass of the host galaxy or 
(in the case of a spiral galaxy) the mass of the stellar bulge. 
Some galaxies, like M32, the dwarf companion to the nearby 
Andromeda galaxy, are known to contain both a steep stel- 
lar cusp (p oc r -1,6 ) and a dynamically-significant black hole 
{M h /M g ~ 0.003). 

The purpose of the present study is to explore the or- 
bital dynamics of realistic triaxial potentials in a system- 
atic way. One goal is to understand how the structure of 
triaxial phase space differs from that of axisymmetric sys- 
tems. Another is to estimate the rate at which chaos implies 
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changes in the structure of a triaxial galaxy, and to un- 
derstand how this rate depends on galaxy properties. In a 
pioneering study, Goodman & Schwarzschild (1981) found 
that the boxlike orbits in a triaxial model with a smooth 
core were often stochastic, but that the behavior of these 
orbits was essentially regular for at least 50 oscillations. On 
the other hand, Merritt & Fridman (1996) found that the 
stochasticity in triaxial models with p cx r~ 2 density cusps 
produced significant changes in the appearance of boxlike 
orbits after just a few tens of oscillations. Merritt & Valluri 
(1996) went on to calculate timescales for mixing in these 
strongly chaotic potentials; they found that ensembles of 
stochastic trajectories could evolve toward invariant distri- 
butions - corresponding to an approximately uniform filling 
of stochastic phase space - on timescale s of only 10 1 — 10 



orbi tal periods. Kandru p and coworkers ( Kandrup & Mahon 
199 4 lM ahon et al- jjgi ) 

noted similar, rapid rates of chaotic 



mixing in a variety of two-dimensional systems. Taken to- 
gether, these results suggest that the characteristic time over 
which chaos manifests itself in the orbital motion of a triaxial 
galaxy is strongly dependent on its radial mass distribution. 
Presumably, this dependence reflects the sensitivity of box- 
like orbits to deflections that occur during close passages to 
the galaxy center. 

Here we use Laskar's frequency mapping technique (§2) 
to look in much more detail at the structure of regular and 
chaotic phase space in realistic triaxial potentials. Laskar's 
technique - summarized in more detail below - allows one to 
quickly and accurately compute the fundamental frequencies 
that characterize the quasi-periodic motion of regular orbits 
on their invariant tori. In addition, one can compute the 
rate at which these frequencies change as a stochastic orbit 
diffuses from one torus to another. The entire phase space 



at a gi 'en energy can then be represented on a single di- 
agram, tne trequency map, wmcn snows tne way m wmcn 



an othdrwise-integrable phase space is modified by the ex- 



istence ot resonances between tne tnree degrees ot freedom. 
The resonances that are most important for determining the 
structure of phase space can immediately be picked out. 

The models considered in this study are the triaxial 
generalizations of the spherical models discussed by Dehnen 
(1993), Carollo (1993) and Tremaine et al. (1994). These 
models have a mass density 



< 7 < 3 



with 



x 2 v 2 z 2 
m =-J + TJ + -J> a>b>c>0, 
cr b z cr 



(1) 



(2) 



and M the total mass. The mass is stratified on ellipsoids 
with axis ratios a : b : c; the x [z] axis is the long [short] axis. 
The parameter 7 determines the slope of the central density 
cusp. For 7 = the model has a finite-density core, while 
for all 7 > the central density is infinite. The strongest 
cusps that we will consider here have 7 = 2, i.e. p cx m -2 
at small radii. Henceforth we will refer to this model as a 
"7-model" or as "Dehnen's model." In the spherical case, 
the radial force in this model is 



9$ 

Or 



GM 



1-7 / r \ 7-3 

1 + - • 

a J \ a/ 



(3) 



which is finite at the center for 7 < 1 ("weak cusp") and 
divergent for 7 > 1 ("strong cusp"). To represent the ef- 
fect on the stellar motions of a central black hole, we also 
consider models with an added central point of mass M%. 
The 7-model is a reasonable description of the stellar den- 
sity at small and intermediate radii in elliptical galaxies. 
Our choice of this model is motivated by the expectation 
that the central regions of a triaxial galaxy are the most im- 
portant for determining the behavior of the boxlike orbits. 
One of our principal objectives will be to establish whether 
there is a critical value of cusp slope 7 or black hole mass 
Mh at which a globally significant transition occurs in the 
structure of phase space. 

We find (§3) that the motion of boxlike orbits in triax- 
ial 7-models is generically chaotic, but that the timescale for 
diffusion in stochastic phase space exhibits great variation, 
both within a single model and as a function of 7 and M%. 
We confirm the result of Papaphilippou & Laskar (1998) 
that the objects of fundamental importance for structur- 
ing phase space are the 3-dimensional resonant tori; in this 
sense, 3DOF systems differ from 2DOF systems, in which pe- 
riodic orbits and their associated families give phase space 
its structure. As the central concentration of the models 
is increased, boxlike phase space becomes more and more 
chaotic. A transition to global stochasticity occurs (§4) when 
the central mass contains between ~ 2% of the galaxy mass; 
for such large values of Mh, the boxlike phase space is al- 
most completely stochastic, and diffusion takes place in a 
very short time. Interestingly, the critical black hole mass 
that we find for transition to global stochastic ity is close to 
the maximum m ass observed in real galaxies (Kormendy & 



Richstone 1995), and to the mass that induces a sud den evo- 



lution toward ax isymmetry in JV-body simulations ( Merritt 
& Quinlan 199£). Central density cusps without black holes 



are less effective at generating chaos; even a modest black 
hole, with Mh ~ 0.3%M g , induces about as much stochastic 
diffusion as the steepest (p cx r~ 2 ) cusps. 

Following Laskar (1993), we characterize the diffusion 
in phase space via the rate of change of the characteristic 
frequencies computed over a fixed interval of time. We find 
that the spectrum of diffusion rates at a given energy can 
be well approximated as a power law with slope close to 
— 1; this spectrum extends over at least six decades (§5). 
Thus, every triaxial potential examined here has a signifi- 
cant population of slowly-evolving stochastic orbits. After 
a long enough time, even weakly stochastic orbits would 
be expected to contribute to chaotic mixing, an essentially 
irreversible process that leads to steady-state, invariant dis- 
tributions in stochastic phase space. Once mixing sets in, 
triaxiality is likely to disappear, since regular box orbits are 
necessary for maintaining triaxial shapes. We estimate (§6) 
the timescale for this to take place as a function of galaxy lu- 
minosity; because typical orbital periods are shorter in faint 
ellipticals, they should be dynamically "more evolved" than 
bright ellipticals and hence less likely to maintain triaxial 
shapes. We predict evolution times that are roughly equal 
to galaxy lifetimes for elliptical galaxies with Mb ~ — 19 to 
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—20, close to the absolute magnitude at which many ellipti- 
cal galaxy properties - including the distribution of intrinsic 
shapes - are observed to change. We are therefore led to a 
picture (§7) in which many of the systematic differences be- 
tween bright and faint ellipticals are a consequence of their 
different dynamical ages. 



with the same axis ratios as the density - that divide the 
model into 21 sections of equal mass. Thus shell 1 encloses 
1/21 of the total mass, shell 2 encloses 2/21, etc.; shell 21 lies 
at infinity. The period of the a;-axis orbit will often be used 
as a reference time; (jjq is defined as 2tt/T x , the frequency of 
the axial orbit. 



2 NUMERICAL TECHNIQUES 

Laskar (1990) developed a new technique, the numerical 
analysis of fundamental frequencies (NAFF), for analyzing 
chaotic motion in the solar system. Laskar's technique is 
based on the principle that an orbit in a system for which the 
motion is close to integrable is described by a well-defined 
set of frequencies uji, one per degree of freedom. For a regular 
orbit, the fundamental frequencies are invariant with respect 
to the coordinates in which the motion is represented. An 
analysis of the variation of the fundamental frequencies with 
respect to time enables one to detect the presence of chaos 
on shorter timescales than are usually possible with more 
classical techniques like Liapunov exponents. In addition, 
the NAFF algorithm allows one to quantify the degree to 
which a stochastic orbit has evolved over a given period of 
time; by contrast, Liapunov exponents reveal only the time- 
averaged rate of divergence in the vicinity of the orbit. A 
plot of the fundamental frequencies of a set of orbits drawn 
from a regular grid of initial conditions yields the "frequency 
map," effectively a representation of the Arnold web of the 
system. The frequency map permits an identification of the 
resonant layers associated with stable and unstable motion 
as well as a determination of the size of the chaotic zones. 

Laskar and collaborators have applied the NAFF al- 
gorithm to a wide variety of systems of astronomical and 
non-astronomi cal interest (|Laskar 199C jLaskar i F^o^schle& 
Celletti 1992t |Laskar 1993; [Papaphilippou fc Laskar 19961 



1998). These authors have shown that the frequency anal- 
ysis technique is a remarkably powerful way to probe the 
structure of phase space in complex, non-integrable systems. 

In the following subsections, we give a brief overview 
of Hamiltonian systems; summarize the basic principles of 
Laskar's technique; and discuss several factors that affect 
the accuracy with which the fundamental frequencies can 
be recovered. 

The models that we consider, both here and below, are 
"maximally triaxial" 7-models, i.e. they satisfy 



(4) 



Henceforth we adopt units in which the total mass M, the 
x-axis scale length o, and the gravitational constant G are 
unity. The mass of the central black hole, when present, 
will be expressed as Mh\ since the total mass of the triaxial 
model is equal to one in our units, Mh may be interpreted as 
the mass of the black hole in units of the galaxy mass. Ex- 
pressions for the gravitational potential and forces are given 
by Merritt & Fridman (1996). Following those authors, we 
assign orbits one of a set of 20 energies, defined as the values 
of the potential on the :r-axis of a set of ellipsoidal shells - 



2.1 Frequencies in Hamiltonian Systems 

If a Hamiltonian system with A degrees of freedom (DOF) is 
integrable, the Hamiltonian H (J, 6*) can be written purely in 
terms of the N actions Jj, H(J,9) = H(Jj). The equations 
of motion of the system are then 

^=0, 0j = §^=Wi(J), 



(5) 



for j = 1,2,..., N. Orbits in the system can be written in 
terms of the complex variables 



(6) 



If du>j(J)/dJj 7^ 0, the system is non-degenerate and in 
principle the actions can be written as 



Jj = Fj(LU ly L02, —, wjv). 



(7) 



When a system is integrable, the motion of a particle pro- 
jected onto the (Jj,0j) plane describes pure circles. This 
is because the motion in phase space occurs on the sur- 
face of tori that are the products of true circles of radii 
\Jj\. The motion around a torus occurs at a rate deter- 
mined by a frequency vector (o>i, UJ2, — , Wjv) which is fixed 
for each torus. In general, we do not know the action-angle 
variables (Jj,9j), but in restricted cases close approxima- 
tions (J'j,8'j) to these variables may be known. In these 
coordinates (J'j,9'j) the motion is no longer pure circles 
but it is still determined by the same vector of frequen- 
cies 9'j = dH/d.J'j = Wj(J) provided the new variables are 
canonically conjugate. The new coordinates may be written 
as a Laurent series expansion in terms of the action-angle 
variables Zj(t): 



C(t) = Zj(t) + 



i(m-w)t 



(8) 



In the limit that the coordinates (JLO'A are action-angle 
variables, the amplitudes a m tend to zero. 

Realistic systems with more than one DOF are rarely 
integrable, i.e. are rarely characterized by N global invari- 
ants. In some cases, the Hamiltonian can be written as a 
perturbation of an integrable Hamiltonian Ho, 



H(3,9)=H (J) + eH 1 (3, 



(9) 



If the perturbation e is small, the KAM theorem guarantees 
that a large fraction of the tori still persist and the motion 
of most of the orbits is still quasi-periodic. The motion on 
this discontinuous set of tori (defined by a Cantor set of fre- 
quency vectors Q = (a>i, u>i, ...,uj n ) and therefore referred to 
as Cantori) is still describable in terms of action-angle vari- 
ables. However it is no longer possible to define a single set of 
action variables that are valid over the entire phase space as 
in equation (7) , since the frequencies change discontinuously 
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wherever the tori are destroyed. Nevertheless, fundamental 
frequencies continue to exist for all regular orbits. 

A central role is played by the resonant tori, tori on 



which he fundamental frequencies satisfy a relation like 
m ■ uj = where m is a vector with integer components. 
The KAM theorem guarantees that "very nonresonant" tori 
- tori for which m • u) is large - persist under small perturba- 
tion of any integrable Hamiltonian. Most of phase space is 
occupied by very nonresonant tori, but resonant tori are also 
dense in the phase space and do not survive small pertur- 
bations. Stable resonant tori (or "elliptic manifolds") gen- 
erate regions of regular motion while unstable resonant tori 
("hyperbolic manifolds") are often associated with stochas- 
tic layers. According to the Poincare-Birkhoff theorem, el- 
liptic and hyperbolic resonances come in pairs, leading to 
alternating regions of regular and stochastic motion. While 
resonant tori are dense in the phase space, only those of suf- 
ficiently low order - i.e., those with |m| sufficiently small - 
are likely to affect more than a very small phase space re- 
gion. Resonances appear as lines in the frequency map, and 
the strength of a particular resonance is immediately appar- 
ent from the degree to wh ich the othe rwise-regular grid of 
points has been distorted (Laskar 1992 ) . 

In a 2 DOF system, motion on resonant tori is closed, 
since the resonance condition implies uj\/u2 = —m/l guar- 
anteeing closure after m revolutions in 9\ and / revolutions 
in 6*2- It is therefore appropriate to state that the phase 
space of a 2 DOF system is structured by its periodic or- 
bits. In a 3 DOF system, however, motion on resonant tori 
is not generally closed, since the single condition m • u = 
does not guarantee that any two of the three frequencies 
are commensurate. The resonance condition guarantees only 
that the trajectory inhabits a submanifold of dimensional- 
ity two on the 3-torus; closure requires the existence of an 
additional, independent resonance relation. While resonant 
tori continue to lend phase space its structure, one can not 
assign every regular orbit in a 3 DOF system to a family 
associated with a particular periodic orbit. 



2.2 Numerical Recovery of the Fundamental 
Frequencies 

Recovering the fundamental frequencies of a regular orbit 
is an aspect of "torus construction," the derivation of the 
map that relates Cartesian and action-angle variables. Two 
general approaches to torus construction h ave been worked 



out in recent years . Iterative techniques (Chapman, Gar- 
rett fc Miller 1976j |Ratcliff, Chang fc Schwarzschild 1984j 



McGill fc Binney 199C|; |Warnock 1991 | [Binncy fc Kumai 
199E ; |Kaasalaiiicn fc Binncy 1994| ) begin by specifying the 



frequencies or actions of an orbit and then attempt to con- 
struct the relations x(#) that must be satisfied if the Oi are 
to increase linearly with time. These methods are often ele- 
gant but suffer from a lack of robustness: the equations to be 
solved are nonlinear and convergence can depend sensitively 
on the accuracy of the initial guess. In a complex dynamical 
system containing many families of orbits, like the triaxial 
potentials considered here, it is generally necessary to "fine- 



tune" the algorithms in order to achieve convergenc e in each 
distinct phase-space region (e.g. Kaasalai ncn 19951) 



Trajectory- follo wing algorithms (Boozer 1982 



Kuo- 



Petravic et al. 1983 ) avoid these problems by making explicit 



use of the quasi-periodic nature of regular orbits. Quasi- 
periodicity guarantees that the motion in any canonical co- 
ordinate x can be expressed as 



x(t) = 2^ a k el 



(10) 



where the uik are linear combinations of the fundamental 
frequencies, u) k = l^>i + muj2 + Wjjs, and the a k are com- 
plex amplitudes. By Fourier transforming the motion, the 
discrete peaks in the spectra can be identified as well as the 
corresponding amplitudes allowing a direct determination of 
x(9). The actions then follow from Percival's (1974) formula, 
completing the specification of the action-angle variables. 
Binney & Spergel (1982, 1984) pioneered this approach in 
the context of stellar dynamics, using a least squares tech- 
nique to fit a model of the frequency spectra to x(t). 

Laskar's (1990) trajectory-following approach, adopted 
here, achieves much greater accuracy. Laskar begins by re- 
placing x{t) in equation ( p"o| ) by f(t), some complex function 
of the orbit, e. g. f(t) = x(t) +iv x (t), which is similar to the 
function £(i) that appears in equation (^). A close approxi- 
mation to f(t), 



/'(*) 



N 

£ 

k = l 



a k e 



(11) 



is then obtained in the following way. First, the time se- 
ries f(t) is translated to an interval [—T/2, T/2] symmetric 
about the time origin. Next, a discrete Fourier transform is 
taken of / and the locations of the peaks ui' k are identified. 
The position of any peak will be defined to an accuracy of 
~ 1/MT where M is the number of equally-spaced inter- 
vals at which / is recorded. The estimate of the location 
strongest peak is then refined by computing the maximum 
of the function 



(w) = </(*). O 



T/2 



f{ty x{t)dt 



(12) 



T/2 



where x(t) — 1 + cos(27rt/T) is the Hanning window func- 
tion. The integral is approximated by interpolating the 
discretely-sampled /. The Hanning filter broadens the peak 
but greatly reduces the sidelobes, allowing a very precise de- 
termination of u)[. Once the first frequency component has 
been found, its complex amplitude is obtained by project- 
ing e 1 ^ 1 * onto the original function f(t), a[ = (f(t), e 1 " 1 '). 
The first term in the series is then subtracted, fi(t) = 
f(t) — aie"'' 1 and the process repeated on fi(t) to find the 
second frequency lu' 2 - In general, the functions ei = e 2 ^ 1 * 
and e 1 " 2 * will not be orthonormal basis functions, but they 
can be made so by applying a Gram-Schmidt orthogonaliza- 
tion process to obtain e2 of the form e2 = bie 1 ^ 2 * — foei, 
with 6i and 62 constants. e2 is then projected onto fi(t) to 
obtain the corresponding amplitude a' 2 . This process is re- 
peated until the residual function does not signficantly de- 
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crease following subtraction of another term. Typically, an 
accurate estimate of fit) is obtained with about 30 terms. 
We have found that extracting a larger number of terms 
from the function, while possible, does not always improve 
the accuracy with which the function f(t) is estimated. 

The frequency analysis of f(t) yields a set of frequen- 
cies uik which are linear combinations of the fundamental 
frequencies (uii, ui2, ^3). Inferring the fundamental frequen- 
cies from the u)u is an ill-defined problem, however, since any 
linearly-independent combination of the fundamental fre- 
quencies can equally well be interpreted as "fundamental." 
Suppose that a particular spike in the frequency spectrum 
has uik = luJi+mu 2 + nu! 3 . One can define new "fundamental 
frequencies" uj[ as 



n 2 ltu' 1 + n 2 2^2 + «23^3> 
n-nu'-t + 71321^2 + n 33 U)' 3 . 



m 
u 3 

In terms of these new frequencies, uJk becomes 
cu k = (/7111 + mn 2 i + nn :il )u'i + 

(Zni2 + mn 2 2 + nn 3 2)w 2 + 

(ini3 + miU3 + nn 33 )uj 3 
= l cui + m L02 + n uj 3 . 

In practice, one expects that the largest amplitude term in 
each frequency spectrum will correspond to a fundamental 
frequency. However, this is only guaranteed to be true if the 
coordinate in question is close to an angle variable. 

This ambiguity turns out not to be serious in the case 
of boxlike orbits, which can be thought of as perturbed rec- 
tilinear orbits. The fundamental frequencies of a box orbit 
correspond approximately to motion in the x, y and z direc- 
tions, and so a frequency analysis of the motion expressed 
in Cartesian variables yields highest-amplitude terms that 
can almost always be identified with fundamental frequen- 
cies. The same is not true for the tube orbits. For instance, 
the highest- amplitude term in both the x- and y-spectra of a 
short (z-) axis tube orbit is typically associated with the fre- 
quency of rotation around the z axis, uj^. The term of next 
highest amplitude typically has a frequency of uir — uj^,, with 
lur the fundamental frequency associated with motion in the 
radial direction. 

Despite these complications, we have successfully used 
an automated, integer-programming "branch and bound" 
algorithm from the NAG library (H02BBF) to extract the 
fundamental frequencies from the frequency spectra. The en- 
tire set of frequencies cut is first sorted in descending order of 
their real amplitudes. The frequency component of highest 
amplitude is defined to be the first fundamental frequency 
u)\. The entire series is then searched in descending order 
of amplitude to find the second fundamental frequency 0^2, 
defined as the first term in the series which is not a simple in- 
teger multiple of oj\ . The integer programming routine then 
defines as ui 3 the next frequency component that is not a lin- 
ear combination of uii and u)2- In the case of box orbits and 
most tube orbits, the error in identifying the fundamental 
frequencies, err(cjfc) = {uj^ — (ZfcWi + nfcW2 + mkU 3 )\, is less 
than about ICP 5 . Papaphilippou & Laskar (1996) showed 



that polar coordinates were better suited than Carestian co- 
ordinates to obtaining the fundamental frequencies of tube 
orbits in an unambiguous way. However, we found that an 
analysis of the tube orbits in terms of Cartesian coordinates 
was satisfactory. Two of the "fundamental frequencies" of 
a tube orbit returned by the integer programming routine 
could almost always be interpreted as ui^, and ujr — lo^, as 
discussed above; the third frequency corresponded to motion 
in a direction parallel to the symmetry axis of the tube. 

Once the fundamental frequencies are obtained, one can 
define the frequency map, a map from initial condition space 
to the space of fundamental frequencies. Since all orbits from 
a given ensemble were selected to have the same energy, 
the distribution of points in frequency space defines a two- 
dimensional surface; we follow the practice of Laskar and co- 
authors of plotting u>i /uj 3 vs. 102/^3 ■ As discussed above, reg- 
ular regions in the frequency map are often associated with 
resonances between the fundamental frequencies, m ■ uj = 0. 
The highest-amplitude terms in the frequency spectrum of 
an orbit near a stable resonance typically reflect the reso- 
nance condition. For instance, near the x — z "banana" or- 
bit, which has uj x /uj z = 1/2, one typically finds u)i/u) 3 = 0.5 
with very high accuracy. An additional frequency, associ- 
ated with the slow libration around the resonance, is also 
present in the spectrum but is generally small (< 10 -6 x u>i) 
and of such low amplitude that it is not selected by the 
integer-programming routine. We again follow the practice 
of Papaphilippou & Laskar (1996, 1998) of defining the fun- 
damental frequencies in the vicinity of a resonance as the 
frequencies corresponding to the highest-amplitude terms. 
Thus, at least two of the fundamental frequencies of orbits 
in an island around a boxlet will remain in exact proportion 
as one moves across the island; the frequency corresponding 
to slow libration is ignored. This practice results in the set 
of orbits around a stable resonance mapping into a line in 
the frequency map, a very convenient representation. 

Stochastic orbits are not quasi-periodic and hence their 
frequency spectra can not be interpreted in terms of just 
three fundamental frequencies. However, if an orbit is only 
weakly stochastic, the integer programming routine will still 
yield a set of "fundamental frequencies" in terms of which 
the various peaks in the frequency spectrum can be ex- 
pressed with more or less accuracy. As the degree of stochas- 
ticity increases, these "fundamental frequencies" become 
simply the three u>k corresponding to the three strongest 
peaks in the frequency spectrum. So defined, the "funda- 
mental frequencies" of a stochastic orbit will depend on the 
time interval over which the orbit was integrated. By inte- 
grating a stochastic orbit over two successive time intervals 
[0, T], [T, 2T], one can compute the change in these frequen- 
cies, w hich is a mea sure of the rate of diffusion in frequency 
space dLaskar 199^ ). 

In the remainder of the paper we will refer to our al- 
gorithm for carrying out the frequency analysis and inte- 
ger programming as "NAFF," even though it differs slightly 
from Laskar's algorithm. 
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Figure 1. Dependence of Laskar's bootstrap accuracy estimate 
8ui = \u>" —u)'\ on the orbital integration time. The solid lines are 
for a regular orbit, with the time series sampled at two different 
intervals At as indicated. The broken lines are for a stochastic 
orbit. 



Figure 2. Diffusion rate Au obtained by measuring the frequen- 
cies of a regular orbit over two successive time intervals of length 
T, plotted as a function of T, for three different values of the 
tolerance parameter of the orbit integrator. 



2.3 Accuracy of the NAFF Algorithm 

In this section, we discuss some of the factors that affect 
the accuracy with which the fundamental frequencies of an 
orbit can be recovered. 

Orbits were integrated using an explicit Runge-Kutta 
scheme of order 8 based on the method of Prince & Dor- 
mand (1981). The program DOP853 from Hairer, Norsett 
& Wanner (1993) uses local error estimation and adaptive 
step size control based on embedded formulas of orders 5 
and 3 respectively. An important feature of the code is its 
ability to produce "dense output," accurate estimates of 
the dependent variables at time intervals much shorter than 
the integration time step. The program DOP853 produces 
dense output via 7th-order interpolation between the output 
points. The dense output option greatly speeded up the con- 
struction of the time series needed for the frequency analysis 
routine. 

In order to estimate the accuracy of the numerically- 
computed u>i, Laskar (1993) suggested the following boot- 
strap scheme. Call /' it) the numerical approximation of the 
original time series f(t), as in equation (11). Carry out a 
second frequency analysis of f'(t). The result is 

JV 

= 5>* eK ''- (13) 
fc=i 

Then 5a k = \a' k — a' k \ and 8u k — \u>' k ' — cu' k \ are estimates of 
the precision with which a' k and u>' k have been determined. 

Laskar (1996) has shown that frequency analysis with a 
Hanning filter yields estimates of the fundamental frequen- 
cies whose accuracies scale asymptotically as 1/T 4 , with T 
the integration interval. For an ordinary FFT the depen- 
dence is only 1/T. Two additional factors affecting the ac- 
curacy of the recovered frequencies are the sampling interval 



At and the accuracy of the numerical integrator. We esti- 
mated the accuracy of the NAFF routine by applying it to 
two orbits in a 7-model with c/a — 0.5 and 7 = 0.5. Fig- 
ure ji] shows how the error in the fundamental frequencies 
of two different orbits varies with T and At; the time unit 
is the period of the long-axis orbit of the same energy. 5ui 
was defined as \cj" — cj'\ for the frequency component of the 
largest amplitude. The solid lines were obtained for a regu- 
lar orbit, and the broken lines were obtained for a stochastic 
orbit. The open circles correspond to a larger sampling in- 
terval of At ~ 10~ 2 (~ 100 steps per orbital period) and 
the filled circles correspond to a smaller sampling interval of 
At ~ 5 x 10 -3 (~ 200 steps per orbital period). This figure 
illustrates the reduction in 8lo for a regular orbit as the inte- 
gration interval is increased. The smaller sampling interval 
also improves the bootstrap estimate of the accuracy 5ui by 
nearly two orders of magnitude. In the rest of the paper we 
adopt an integration interval of 50 orbital periods (defined 
as the period of the long-axis orbit of the same energy) and 
a sampling interval At = 5 x 10 -3 . Figure [j] also shows that 
the "frequencies" of a stochastic orbit are obtained with a 
much lower accuracy (10~ 4 — 10 -6 ) than for a regular orbit. 
Clearly, the function f'(t) obtained for a stochastic orbit is 
not a good approximation to the original f(t). 

Another way to estimate the accuracy of the NAFF rou- 
tine is to integrate a single, regular orbit over two successive 
time intervals and compare the estimates of the fundamental 
frequencies. We define Auj — \u(Ti) —u){T2)\/u>qT, where Ti 
and Tb are the midpoints of the two intervals, T = T2 — Ti, 
and ojo — 2n/Tx, the frequency of the long-axis orbit of the 
same energy. Figure ^ shows Alo as a function of T and of 
the tolerance parameter of the numerical integrator for the 
regular orbit of Figure |lj. Figure ^| indicates that decreasing 
the tolerance parameter below 1 x 10 -8 does not significantly 
affect the value of Auj when the sampling interval is fixed to 
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this value. Unless otherwise indicated, results in this paper 
were obtained with a tolerance parameter of 10 -8 . 



in this potential ( [Fridman fc Merritt 1997[ ) . Instability of the 
z— axis orbit first occurs at low energy through bifurcation 
of the 1 : 1 y — z loop orbit (Goodman & Schwarzschild 



3 THE STRUCTURE OF PHASE SPACE 

We begin by looking in detail at a triaxial 7-model with 
c/a — 0.5 and 7 = 0.5, a "weak cusp." Two ensembles of ~ 
10 4 orbits were integrated in this model at the fixed energy 
corresponding to shell 8, which lies just inside the half-mass 
radius. Initial conditions for the two ensembles were chosen 
so as to generate either boxlike orbits or tube orbits; i.e. 
orbits with stationary points, or orbits that circulate about 
one of the principal axes of the model. (A more complete 
discussion of the two initial-condition spaces may be found in 
Merritt & Fridman 1996). Orbits were integrated for a total 
of 100 periods of the a;-axis orbit. Fundamental frequencies 
were computed independently for the first and second halves 
of this interval so that diffusion rates could be calculated, 
as described above. 



3.1 Box-Orbit Initial Conditions 

A large part of the phase space of a triaxial system is as- 
sociated with orbits that have a stationary point, i.e. that 
touch the equipotential surface with zero velocity. All of the 
box orbits in an integrable (Stackel) potential fall into this 
category, as do the "centrophobic" boxlets, like the x — z ba- 
nana, in non-integrable potentials. ( "Centrophilic" boxlets, 
like the anti-banana, typically do not have a stationary 



point. Such orbits are generally unstable (Miralda-Escude 
fc Schwarzschild 1989| )). In Figure |, the NAFF algorithm 
has been applied to a set of 9408 orbits whose initial condi- 
tions lay in a regular grid on one octant of the equipotential 
surface corresponding to shell 8. Figure presents a pla- 
nar projection of the initial condition space; the gray scale 
has been adjusted in proportion to the logarithm of the dif- 
fusion rate, defined as the change in the largest-amplitude 
fundamental frequency as evaluated over the initial and fi- 
nal intervals of 50 periods each. Figure ^b is the frequency 
map corresponding to the first integration interval. In both 
figures, resonances defined by lw x + mu y + nui z = have 
been identified for certain integer values of (l,m,n). 

The frequency map exhibits some regions in which the 
points are arranged in a more-or-less orderly fashion. These 
regions may be identified with parts of phase space that 
have retained their regular character in spite of the perturba- 
tion of the potential away from integrable, or Stackel, form. 
Separating these regular regions are the resonance lines: ei- 
ther solid lines, corresponding to stable resonant tori; or 
empty gaps, corresponding to stochastic layers. In addition, 
a number of points are scattered in a more irregular manner 
around the figure. These are the stochastic orbits. However, 
even the stochastic orbits are mostly clustered around res- 
onance lines, indicating that the stochastic motion is often 
strongly influenced by structures in phase space that are 
similar to invariant tori. 

The short (z— ), intermediate (y— ) and long (as—) axis 
orbits are all unstable to lateral perturbations at this energy 



1981). The influence of this resonance is visible at the top 
of Figure ^b, where many of the stochastic orbits started 
near the z axis are clustered around the (0, 1, —1) resonance 
line. (The stable branch of the bifurcation, the y — z loop 
orbits, do not appear in this plot since they have no sta- 
tionary points.) Instability of the y— axis orbit first occurs 
through the 1 : 1 x — y loop bifurcation; the region around 
the (1,-1,0) resonance is also well populated in the fre- 
quency map, by stochastic orbits that start near the y— 
axis. Finally, the x— axis orbit becomes unstable at the bi- 
furcation of the 1 : 2 x — z banana orbit, which occu rs at 
about shell 2 in this model ( |Fridman fc Merritt 1997[ ). Un- 
like the 1 : 1 loop orbits, however, the stable banana orbit 
does lie in stationary initial condition space, and hence the 
family of stable orbits generated by the bifurcation appear 
in the frequency map along the (2, 0, —1) resonance line. A 
smaller number of stochastic orbits lie near to, but offset 
from, the (2,0, —1) line. 

As emphasized by Papaphilippou & Laskar (1998), the 
motion in boxlike phase space is strongly influenced by res- 
onances between the three degrees of freedom. Such reso- 
nances are dense in the phase space of an unperturbed, i.e. 
integrable, Hamiltonian, in the sense that every torus lies 
near to a torus satisfying m ■ u — for some (perhaps very 
large) integer vector m. However, most tori are very non- 
resonant in the sense that m • w is large compared with 
| m |-(jv+i)^ w jth TV the number of degrees of freedom. These 
very non-resonant tori persist under finite perturbations of 
the Hamiltonian and may be identified with the orderly 
regions between the resonance lines in the frequency map 
of Figure ^b. Resonant tori, particularly those with small 
m, are either destroyed by a finite perturbation, generating 
stochastic motion, or become associated with a regular re- 
gion surrounding the resonance, or (typically) both. Many - 
perhaps most - of the regular orbits in Figure tm can be iden- 
tified with a stable resonance; among the most important are 
the (2,1,-2), (3,-1,-1), and (4,-2,-1) resonances. The 
intersection of any two resonance lines defines a periodic or- 
bit, i.e. an orbit for which uj x /lu z = l/n and uj y /uj z = m/n. 
When stable, a periodic orbit can have its own associated 
regular region, distinct from the regular regions associated 
with the 3D resonances. At least two such regions are appar- 
ent in Figure associated with the 5:6:8 and 7 : 9 : 10 
periodic orbits (marked 1 and 2). In addition, the x — z 
fish (2 : 3) and x — y pretzels (3 : 4) planar boxlets have 
associated regular regions. 

It is sometimes stated that all of the regular orbits in a 
galactic potential belong to families that can be associated 
with stable periodic orbits. For instance, the box orbits in 
a Stackel potential are commonly associated with the long- 
axis orbit, while the tube orbits are said to belong to fam- 
ilies generated by the x — y or y — 2 1 : 1 planar loops. 
While the identification of regular orbits with stable peri- 
odic orbits is sometimes justified, only a small fraction of 
the regular orbits in Figure ^|a can be usefully assigned to 
families associated with periodic orbits. A larger fraction of 
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Figure 3. Box-orbit phase space at shell 8 in the model with 7 = 0.5. (a) Initial condition space: one octant of the equipotential surface 
has been projected onto a plane. Each orbit begins on this surface with zero velocity. The top, left and right corners correspond to the 
z (short), x (long) and y (intermediate) axes. The grey scale is proportional to the logarithm of the diffusion rate of orbits in frequency 
space; initial conditions corresponding to regular orbits are white. The regions labelled "1" and "2" contain orbits associated with the 
5:6:8 and 7 : 9 : 10 periodic orbits, respectively. Other regions are labelled with the integers (I, m, n) that define resonance zones, (b) 
Frequency map: the fundamental frequencies returned by the NAFF algorithm are plotted as rotation numbers lu x /u> z and uj y /u z for 
each of the orbits in (a). The most important resonances are labelled. Stable resonances produce solid lines; gaps correspond to unstable 
resonances. 



the regular orbits are located in resonance zones, regions as- 
sociated with tori that satisfy a single resonance condition 
lu>x +rncj y + nuj z = 0. We expect orbits in such regions to be 
approximately confined to two-dimensional sub-manifolds; 
in configuration space, they should take the form of thin, 
curved sheets. An example of such an orbit, which we call a 
"thin box," is shown in Figure ^| Closure requires an addi- 
tional resonance condition m • u = to be satisfied so that 
the frequency vector can be written as u — ncoo with n an 
integer vector. 

In initial condition space (Figure [^) , the regions asso- 
ciated with stable resonance zones appear as bands of white, 
flanked on both sides by narrower stochastic strips. Periodic 
orbits lie at the intersection of two or more of these bands, 
and some of these periodic orbits - e.g. the 5:6:8 and 
7 : 9 : 10 orbits in Figure ^a - are surrounded by regions 
in which the motion is regular. Orbits in these restricted re- 
gions may be usefully associated with the periodic orbit at 
their center. But Figure ^a suggests that the regular regions 
associated with stable periodic orbits are relatively small, 
and they are often separated from the larger regular regions 
by strips where the motion is stochastic. There would seem 
to be no reason to associate the regular motion in the largest 
part of any resonance zone with any particular periodic or- 
bit. 



In a 2 DOF system, the condition Iwi + mu>2 ~ that 
defines a resonant torus also defines a periodic orbit, since 
W1/W2 = —m/l implies commensurability of the motion in 
the two directions around the torus. The distinction between 
periodic orbits and resonant tori therefore does not appear 
in systems with fewer than three degrees of freedom. In the 
same way, it is clear from Figure that orbits restricted to, 
say, the x — y plane are affected by the (4, —3, 0) resonance 
zone only where that zone intersects the x — y plane, at 
the x — y pretzel orbit. In the case of the majority of non- 
planar regular orbits, however, no unique identification with 
a periodic orbit is possible. 

In addition to the regular orbits associated with reso- 
nance zones and with periodic orbits, a third class of regu- 
lar orbit is apparent in Figure |S| These are orbits that he 
in the (essentially) completely regular regions that lie be- 
tween the resonance zones. Orbits in these regions may be 
identified with very non-resonant tori, i.e. tori that have 
maintained their integrity in spite the perturbation of the 
Hamiltonian away from integrable form. In an integrable po- 
tential, all of the regular orbits would belong to such regions; 
in a non-integrable potential, only tori that are sufficiently 
non-resonant are expected to avoid being strongly deformed 
by the perturbation. Strictly speaking, we expect to find 
some high-order resonance zones and associated stochastic 
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Figure 4. A thin box orbit associated with the (2, 1, —2) res- 
onance. The lower view is a cutaway showing that the orbit is 
confined to a membrane. 



regions imbedded within even these "regular" regions, and in 
fact a close inspection of Figure reveals just such struc- 
ture. However, the chaotic diffusion rates in these regions 
are generally small and for practical purposes the motion in 
them is effectively regular throughout. Again, we note that 
the motion in one of these regular regions can not be prop- 
erly identified with any single periodic orbit; instead, the 
regularity of the motion would seem to be properly identi- 
fied with the persistence of integrability in extended, non- 
resonant parts of phase space in spite of the perturbation of 
the Hamiltonian away from fully integrable form. 

An equivalent way to define the three "types" of regular 
orbit described above is in terms of the number of indepen- 
dent relations m • u) = that characterize their associated 
phase space regions, i.e. in terms of the degeneracy of the 
resonance. A periodic orbit is defined by two such relations; 
a resonance zone by one; and the regular regions that persist 
in spite of the perturbation of the potential away from inte- 
grability are defined by none. It is clear from this definition 
that the number of distinct "types" of regular orbit is equal 
to the number of degrees of freedom. While regular orbits 
of all three types occupy three-dimensional regions in con- 
figuration space, orbits associated with a single resonance 
condition are approximately confined to sheets, while orbits 
associated with two resonance conditions are approximately 
confined to closed curves. 

Given that essentially none of the boxlike orbits in Fig- 



ure ga can be properly identified with the long-axis orbit, 
it is reasonable to ask in what sense the box orbits in any 
triaxial potential - even that of the Perfect Ellipsoid - are 
correctly said to belong to the "family of orbits" associated 
with the long-axis orbit. We would argue instead that the 
regularity of the box orbits in a Stackel potential is prop- 
erly ascribed to the global integrability of the potential, and 
not to the local stability of any single periodic orbit like the 
long-axis orbit. 

The regions corresponding to stochastic motion in Fig- 
ure can likewise be separated into three fairly distinct 
types, with one important qualification. Stochastic orbits 
live in a phase space of higher dimensionality than regu- 
lar orbits, and in fact we expect every stochastic orbit to 
eventually visit (via Arnold diffusion) every point of phase 
space, even points that are arbitrarily close to regular or- 
bits. However the timescale for Arnold diffusion is extremely 
long, and for practical purposes, stochastic orbits - like reg- 
ular orbits - can often be associated with restricted parts 
of phase space. Accordingly, in Figure ^a, we identify three 
"types" of stochastic orbit. First are the stochastic orbits 
closely associated with an unstable periodic orbit. An ex- 
ample is the region around the 4:5:6 periodic orbit that 
lies at the intersection of the (3,0, —2) and (4, —2, —I) res- 
onance zones. Second, some stochastic orbits lie along the 
edges of the stable resonance zones, as mentioned above. 
Third, there is a large stochastic strip connecting the y— 
and z— axis orbits. Goodman & Schwarzschild (I98f) first 
described this stochastic strip in their study of stochastic 
motion in a triaxial model with a core. The overlap of the 
chaotic zones associated with the y and z-axis orbits appears 
to be responsible for most of this strip, although a number 
of other resonances are also important, including especially 
the (f, —2, 1) resonance. Although the diffusion rate of or- 
bits in the stochastic strip is high, we would not expect these 
orbits to wander ergodically over the entire energy surface 
due to the existence of substantial regular regions which 
would inhibit the diffusion. However, the near constancy of 
the diffusion rate throughout the y — z strip suggests that 
the motion is nearly ergodic over this restricted region, at 
least for elapsed times of ~ I0 2 oscillations or more. 



3.2 Tube-Orbit Initial Conditions 

The other major category of orbits in a triaxial potential, 
the tube orbits, are characterized by a finite, time-averaged 
angular momentum about the long or short axis. In an in- 
tegrable triaxial potential, the three families of tube orbits 
- the short-axis tubes, the inner long-axis tubes, and the 
outer, long-axis tubes - can all be recovered by taking ini- 
tial conditions in the x — z plane, with v y = ( [gchwarzschild 



1992 ) . This initial condition space will also include some box 
orbits. Figure [B^ presents the x — z initial condition space 
for the same model of Figure ^; again, the grey scale is ad- 
justed in proportion to the logarithm of the diffusion rate. 
The corresponding frequency map is shown in Figure [sja. 
Any non-periodic orbit intersects the x — z plane in at least 
two points; in the case of tube orbits, the two points lie in- 
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Figure 5. Tube-orbit phase space at shell 8 in the model with 7 = 0.5. (a) Initial condition space: orbits begin on the x — z plane 
with v x = v z = 0. The grey scale is proportional to the logarithm of the diffusion rate of orbits in frequency space; initial conditions 
corresponding to regular orbits are white. The regions labelled "S," "O," "I" and "B" contain short-axis tubes, outer long-axis tubes, 
inner long-axis tubes and boxes, respectively. The three most important resonance zones are labelled; they intersect at the (unstable) 
4:5:6 orbit (Figure ^j). (b) Frequency map, with the important resonances labelled. Only the tube orbits are included, lu^ is the 
fundamental frequency associated with motion about the short axis (short-axis tubes), or long axis (long-axis tubes). 



side and outside the curve that defines the corresponding 
"thin" orbit family. 

As discussed above, the highest-amplitude components 
in the x— and y— spectra of a short (z— ) axis tube are typ- 
ically both associated with the frequency of rotation about 
the z axis, uj^\ the frequency corresponding to radial oscil- 
lations, uir, appears typically in terms of much lower ampli- 
tude. For a few of the ~ 10 4 tube orbits integrated here, the 
automated routine for extracting fundamental frequencies 
appeared to fail, leaving some random gaps in the frequency 
map of Figure flb. (Papaphilippou & Laskar (1996) reported 
similar difficulties in their analysis of tube orbits.) However, 
the fundamental frequencies of the great majority of regular 
tube orbits were successfully recovered. In Figure |H|b, only 
orbits that clearly belonged to one of the three tube fami- 
lies are plotted in the frequency map; the boxlike orbits are 
omitted for clarity. The fundamental frequency labeled u;</> 
is defined as the frequency associated with motion around 
the z axis in the case of the short-axis tubes (labeled "S" 
in the figure), and as the frequency associated with motion 
around the x axis in the case of the inner (I) and outer (O) , 
long-axis tubes. 

This part of phase space is largely regular. Aside from 
the box orbits, which are confined to regions near the zero- 
velocity curve, all stochastic tube orbits are associated with 
a handful of narrow resonance zones. These zones exist at 



the transition points between the various families, as noted 
by Merritt & Fridman (1996) and Papaphilippou & Laskar 
(1998). A few additional resonances are important, includ- 
ing the (4, —2, —1), (3, 0, —2) and the surprisingly high-order 
(15, 12, —20) resonance. These three resonances intersect at 
the 4:5:6 periodic orbit, an unstable tubelet (Figure ||). 
The diffusion rates appear to drop to zero along the curve 
defining the thin tube orbits belonging to the S and O fam- 
ilies, a reasonable result. 

We have encountered little discussion in the literature of 
stochastic tube orbits; the nearest examples that we know of 
are orbits that Kandrup, Eckstein & Bradley (1997) describe 
as appearing "alternately boxy and loopy." Accordingly, we 
illustrate in Figure the evolution of a tube orbit associ- 
ated with the unstable (4, —2, —1) resonance zone in Figure 
|B^,. The orbit was integrated for three successive intervals 
equal to 50 periods of the long- axis orbit. The stochasticity 
is evident in the gradual change of the orbit's shape. Assum- 
ing that the entire stochastic phase space at this energy is 
interconnected, we would expect this tube orbit to eventu- 
ally find its way into boxlike phase space. However the time 
required is likely to be extremely long, since the diffusion 
would have to occur along the (4,-2,-1) resonance line, 
and diffusion along resonance layers - i.e. Arnold diffusion 
- is known to be extremely slow. 
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7 close to zero should be close to integrable (though never 
exactly so, since the 7-models do not reduce to integrable 
form even for 7 = 0), while values of 7 greater than about 



one are known to gene rate strongly stochastic motion (Mer- 
ritt & Fridman 1996). The mass of a central singularity. 



Figure 6. The unstable, 4 : 5 
appears in Figure ^ 



resonant tube orbit that 




and the energy, should also act as perturbation parameters; 
for instance, the long-axis orbit generally becomes unstable 
at sufficiently high energies even in triaxial potentials with 
cores. In terms of any of these perturbation parameters, box- 
like phase space should be most strongly affected since tube 
orbits avoid the destabilizing center. 

In the case of extreme perturbations, dynamical sys- 
tems sometimes exhibit a transition to global stochasticity, 
a regime in which the stochastic motion is interconnected 
over large portions of the phase space. (Strictly speaking, 
stochastic phase space in a 3 DOF system is always inter- 
connected through the Arnold web, but Arnold diffusion is 
extremely slow unless the stochastic regions overlap.) In the 
globally stochastic regime, there are few barriers to the mo- 
tion and orbits can wander over the energy hypersurface 
in little more than an orbital period. One standard defini- 
tion of the transition to global stochasticity is based on the 
overlap of t he stochastic la yers associated with the primary 
resonances (Chirikov 196C| ). Other definitions (Lichtenberg 
fc Lieberman 1992[ ) are based on the fraction of phase space 



that is stochastic, the degree to which the stochastic regions 
are interconnected, etc. 

In triaxial stellar systems, a transition to global stochas- 
ticity implies the replacement of distinguishable box orbits 
by orbits that move nearly ergodically over the energy sur- 
face and that densely fill the configuration space volume 
defined by an equipotential surface. Since box orbits with 
various shapes are always found to be strongly popu l ated in 
self-c onsistent triaxial models (Schwarzschild 197£; Statler 
1987), significant triaxiality should be difficult or impossi- 
ble to maintain in a system where the motion is globally 
stochastic. This hypothesis has re ceived some support from 
self-consistent equilibrium studies (Schwarzschild 1993 ; |Mer- 
ritt 1997 ) and from iV-body integrations ( |MerTrtT^T^uirilan 



Figure 7. A stochastic, short-axis tube orbit. The initial con- 
ditions lie in Figure |j]a on the resonance zone labelled (4, —2, 1), 
slightly above the 4:5:6 periodic orbit, (a), (b) and (c) show 
three successive integrations of the orbit, each over a time span 
of 50 periods of the long-axis orbit. 



4 THE TRANSITION TO GLOBAL 
STOCHASTICITY 

When an integrable potential is perturbed by increasingly 
large amounts, the fraction of phase space associated with 
stochastic motion and the rate of stochastic diffusion typi- 
cally increase. In triaxial 7— models, the parameter 7 that 
defines the steepness of the central density cusp is expected 
to play the role of a perturbation parameter: models with 



We investigated the approach to global stochasticity in 
the triaxial 7-models by applying the NAFF algorithm to 
a set of models with different values of the central density 
slope 7 and with central point masses Mh of various sizes. 
(Dependence of the degree of stochasticity on the energy 
is discussed in the following section.) In each model, we 
integrated ensembles of ~ 10 4 orbits at shell 8 in boxlike 
initial condition space and computed their fundamental fre- 
quencies over two adjacent time intervals with the NAFF 
algorithm. We carried out similar experiments in the initial 
condition space of tube orbits; however, as expected, the 
degree of stochasticity was not found to depend strongly on 
the central density structure and so those results will not be 
presented here. 

Figures ^ and ^ show the diffusion rates and frequency 
maps for ensembles of orbits with boxlike initial conditions 
in a set of 7-models with five different values of 7 between 
and 2 and c/a — 0.5. The density cusps in real elliptical 
galaxies are almost always found to have slopes that lie in 
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Figure 8. Initial condition space at shell 8 in a set of models with various values of the cusp slope 7. The grey scale is proportional to 
the logarithm of the diffusion rate in frequency space, defined as the fractional change in the largest-amplitude fundamental frequency 
over 50 periods of the long-axis orbit. White regions correspond to regular orbits. 
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Figure 9. Frequency maps for the orbits of Figure ^. 
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Figure 10. Initial condition space at shell 8 in a set of models with c/a = T = 7 = 0.5, and various values of M^, the mass of a central 
singularity. The grey scale is proportional to the logarithm of the diffusion rate in frequency space, defined as the fractional change in 
the largest-amplitude fundamental frequency over 50 periods of the long-axis orbit (in the model without a central point mass). White 
regions correspond to regular orbits. 
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Figure 11. Frequency maps for the orbits of Figure llOl 
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this range (Merritt fc Fridman 1995; Gebhardt et al. 1996) 



with the steeper cusps usually appearing in fainter galaxies. 
The diffusion rates have been normalized by dividing the 
change in uj by the frequency of the long-axis orbit in all of 
the models. Since the orbital integrations were carried out 
for a fixed number of long-axis orbital periods, the diffusion 
rates plotted in Figure |§] may be interpreted as the fractional 
change in ui that takes place over a fixed number (50) of 
orbital oscillations. 

For 7 = 0, both the diffusion plots and the frequency 
maps reveal that the motion is largely regular, at least over 
the interval of ~ 10 2 oscillations for which the orbits were 
integrated. The stochastic orbits are mostly identified with 
the instability strip that connects the y and z axes in ini- 
tial condition space. A number of orbits that lie close in 
frequency space to a resonance between the three degrees 
of freedom exhibit weaker stochasticity; the most prominent 
resonance is the (1, —2, 1), which extends from the y — z in- 
stability strip through the unstable, 4:5:6 and 3:4:5 
periodic orbits. However the diffusion rates for most of these 
orbits are lower than that in the y — z instability strip by 
orders of magnitude. As a result, the majority of orbits out- 
side of the y — z instability strip remain in a smooth grid in 
the frequency map. 

As 7 is increased from to 1, the area in initial condi- 
tion space associated with resonance layers becomes larger. 
Almost all of the regular orbits in the 7=1 model can be 
associated with a stable resonance, and a smaller number 
with stable periodic orbits that lie at the intersection of two 
or more resonances. The regions between the resonance lines 
in the frequency map now show almost no hint of a regular 
grid; in other words, essentially none of the phase space in 
the 7 = 1 model can be usefully identified with the inte- 
grable phase space of a Stackel potential. Nevertheless, the 
majority of the strongly stochastic orbits remain confined to 
the y — z instability strip; diffusion rates in other parts of 
initial condition space are generally much lower. 

The most dramatic change in the structure of phase 
space occurs as 7 is increased from 1 to 2. The y — z in- 
stability strip enlarges to include most of the equipotential 
surface; the diffusion rates in this interconnected stochastic 
region become nearly constant, implying that the stochastic 
motion is essentially free to explore the entire phase-space 
region accessible to it. Most of the stable resonance zones 
disappear, leaving finally only the (2, 0, —1) and (3, — 1, —1) 
resonances with appreciable numbers of associated regular 
orbits. In fact a large fraction of the regular orbits in the 
7 = 2 model appear to be closely associated with just three, 
stable periodic orbits: the 1 : 2 x — z banana orbit (the 
2,0,-1 resonance); and the 4:5:7 and 5:7:8 orbits 
(both associated with the 3,-1,-1 resonance). Essentially 
all of the points on the frequency map that lie away from 
these two resonance lines are stochastic. However, a certain 
amount of structure remains in stochastic phase space, since 
many of the fundamental frequencies returned by the NAFF 
algorithm for the stochastic orbits are clustered near one or 
more resonance lines. In other words, over time invervals of 
~ 10 2 oscillations, many of the stochastic orbits continue to 



behave to a certain extent like regular orbits, with more-or- 
less well-defined frequencies. 

Figures ^ and [H] present the analogous results for en- 
sembles of orbits in models with 7 = 0.5 and with five differ- 
ent values of Mh, the mass of a central singularity in units 
of the total mass of the model. The progression from es- 
sentially regular to essentially chaotic motion is now more 
dramatic. Even for Mh = 0.0003, the regular grid in the 
frequency map is destroyed, and the motion is dominated 
by the resonances. As Mh approaches 0.01, the y — z in- 
stability strip enlarges to include most of the equipotential 
surface, and the diffusion rates are uniformly high through- 
out most of stochastic initial condition space. The structure 
in the frequency map is almost completely destroyed when 
Mh = 0.03. The only remaining regular orbits are associated 
with the x — z banana orbit and the 3:4:5 periodic orbit, 
while the stochastic orbits are distributed almost randomly 
around the frequency map. In terms of the definitions pre- 
sented above, the phase space of the Mh = 0.03 model (at 
least, that part of phase space associated with boxlike or- 
bits) appears to be in the globally stochastic regime. While 
the precise value of Mh corresponding to the transition to 
global stochasticity is not clearly defined, Figures |Io| and |ll| 
suggest that it lies between 0.01 and 0.03 at this energy. 

By comparison, even the steepest cusp (7 = 2) appears 
to generate no more stochasticity than a central singularity 
with Mh ~ 0.003. In other words, in models without a cen- 
tral black hole, the transition to global stochasticity appears 
to be just beginning as 7 increases above 2. 

The transition to global stochasticity seen in the box- 
like orbits is in sharp contrast to the behavior of the tubes. 
Since the time-averaged angular momenta of tube orbits is 
nonzero, tube orbits never go close enough to the center to 
experience the effects of a steeply-rising central force. Conse- 
quently, one does not expect the structure of phase space as 
seen in either diffusion rate plots or the frequency map (Fig- 
ure^) of tube orbits to show much dependence on cusp slope. 
This was in fact found to be the case. The stochasticity in 
the tube orbits arises from the inherent non-integrability of 
the model, and its character is determined by the overall 
shape of the potential, not by any peculiarity of the central 
density structure. 

It was remarked above that regular orbits could be 
classified in terms of their degree of degeneracy, i.e. in 
terms of the number of independent resonance conditions 
Zu;i + mu)2 + nus = that define their associated phase- 
space regions. Figures ^ and [l(] suggest that a systematic 
change takes place in the typical degree of degeneracy of 
the regular orbits as the strength of the perturbation is in- 
creased. In weakly chaotic potentials (e.g. 7 = 0, Mh = 0), 
most of the regular orbits lie in regions of zero-fold degen- 
eracy, i.e. strongly non-resonant regions that have retained 
their integrable character in spite of the perturbation of the 
Hamiltonian away from integrable form. As the perturbation 
is increased (7 = 0.5), a larger fraction of the regular orbits 
lie in resonance zones, regions defined by a single resonance 
condition. When the perturbation is large (7 = 2), the only 
remaining regular orbits lie in regions of two-fold degener- 
acy, i.e. regions associated with stable periodic orbits. 
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Figure 12. Spectra of diffusion rates for two ensembles of boxlike 
orbits at shell 8 in the model with 7 = 1. The thin curve was 
derived from ~ 10 4 orbits integrated for 100 orbital periods; the 
thick curve is from ~ 10 3 orbits integrated for 400 periods. Tick 
marks separate regular from stochastic orbits, as discussed in the 
text. The thin straight line has a logarithmic slope of —1. 



5 DIFFUSION 

5.1 Spectra of Diffusion Rates 

While the distinction between regular and stochastic motion 
is unambiguous in principle, stoc hastic orbits can m imic reg- 
ular orbits for very lon g ti mes (Contopoulos 1971). This is 
shown clearly in Figure [13, a plot of the distribution of dif- 



fusion rates for boxlike orbits at shell 8 in the 7 = 1 model. 
The thin curve was derived from the same ~ 10 4 orbits 
whose frequency map is displayed in Figure []. The solid 
curve is from a smaller set of ~ 10 3 orbits integrated for a 
total of 400 long-axis orbital periods, four times the inte- 
gration interval of the larger ensemble. (The accuracy of the 
numerical integration routine was increased for the longer 
run.) Since the precision with which the NAFF algorithm 
computes the fundamenta l frequencies is a strong function 
of the integration interval (Laskar 1996), we expect the solid 
curve to extend to smaller values of Auj than the thin curve, 
assuming of course that very weakly chaotic orbits exist in 
this potential. The abscissa, lAw/wolso, is the change in fun- 
damental frequencies over an interval of 50 long-axis orbital 
periods, normalized to ujq, the frequency of the long-axis or- 
bit at shell 8. In the longer run, the measured changes in uj 
were scaled to the shorter interval following the procedure 
described in §6. 

Th e spectrum of diffusion rates is well described as a 



denned break in the spectrum at low Auj. Presumably, many 
of the orbits with Auj below the minimum value detectable 
by the NAFF routine are regular. Based on the accuracy 
tests described above, we estimated, for the shorter run, 
the smallest significant value of Aw. This value, indicated 
by the thin vertical tick mark, implies that ~ 22% of the 
orbits in the shorter run are effectively regular. A similar 
tick mark has been placed on the solid curve at the 22% 
point. Happily, on both curves, the tick mark falls roughly 
at the point where the distribution turns upward: sharply 
in the case of the longer run, more weakly in the case of the 
shorter run. Plots of the starting points on the equipotential 
surface of orbits to the left of the tick marks also show a very 
similar distribution for the two ensembles. We conclude that 
roughly 1/5 of the orbits in each ensemble are regular, or so 
weakly stochastic that we can not distinguish them from 
regular orbits. 

One might reasonably ask whether there are any regular 
orbits in these ensembles. Perhaps the upturn at small Aui 
in Figure |l^ reflects a population of orbits with extremely 
low but nonzero diffusion rates. We consider such an inter- 
pretation to be unlikely. Certain of the periodic orbits in this 
initial condition space (e.g. th e x—z banana) are known to be 
stable to small perturbations (Fridman & Merritt 1997) and 



must have families of associated regular orbits. The orbits 
with the lowest diffusion rates do in fact have initial condi- 
tions lying close to these stable periodic orbits. A number 
of experiments reassured us that the upturn at small Au; in 
the distribution of diffusion rates always occurs at roughly 
the same percentile point in a given initial-condition space, 
regardless of (modest) changes in the integration interval or 
accuracy of the numerical integrator. If there were no bona- 
fide regular orbits, we would expect to see the fraction of 
orbits with significantly nonzero Aw's increase monotoni- 
cally with increasing accuracy of the numerical routines. 

In any case, we clearly see a smaller population of reg- 
ular box orbits here than was seen in many earlier studies. 
For instance, Merritt & Fridman (1996) used Liapunov expo- 
nents to study orbits in the same 7 = 1 model investigated 
here. In their ensemble of 192 boxlike orbits from shell 8, 
they estimated that ~ 40% were regular (their Figure 7a), 
roughly twice the fraction of regular orbits found here. This 
difference is not surprising given the greater sensitivity of 
the NAFF algorithm to weak stochastiticy. Figure jl2] im- 
plies that Merritt & Fridman (1996) were unable to detect 
stochasticity in orbits with Auj/ujq less than about 3 x 10 -4 . 
As we show below, such small amounts of diffusion corre- 
spond to very nearly regular behavior. 

The existence of a population of stochastic orbits with 
very low diffusion rates is consistent with a number of stud- 
ies of motion in "decomposable" phase spaces, i.e. sys tems 
containin g both regular and stochastic trajectories (e.g. Kar- 



power l ^iw in Auj with logarithmic slope close to —1. In the ney 1983] ). In decomposable systems, stochastic orbits are 



run with the longer integration interval (solid curve), this 
power law distribution extends over at least six decades in 
Aui, a remarkable range. Evidently, an appreciable fraction 
of the stochastic orbits in this model exhibit motion that 
is very nearly, though not quite, regular. This fact makes it 
difficult to decide which orbits are regular - there is no well- 



hampered in their diffusion by invariant tori; when a large 
fraction of phase space is associated with such tori, many 
of the stochastic orbits exhibit quasi-regular behavior over 
long periods of time. 

The approximately 1 / Auj distribution of diffusion rates 
in Figure Il3 can not extend to very large or very small Aw's 
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Figure 13. Spectra of diffusion rates for the ten ensembles of 
boxlike orbits, (a) Mh = 0; 7 = 0, 0.5, 1, 1.5, 2, increasing upward, 
(b) 7 = 0.5; M h = 0.0003,0.001,0.003,0.01,0.03, increasing up- 
ward. The thin straight lines have power-law slopes of —1. The 
curves are offset vertically with respect to each other by one unit 
in the ordinate. 

without producing a logarithmic divergence in the total 
number of orbits. In the spectrum derived from the shorter 
run, we in fact observe a dropoff for |Aw/u;o|50 > 0.05. 
This is expected, since the fundamental frequencies asso- 
ciated with orbits at a given energy are restricted to a range 
of values near the oscillation frequencies along the major 
axes, which means that Aui/loq can never exceed unity and 
will generally be much smaller. There must also be a falloff 
for small diffusion rates, but this apparently occurs at such 
small values of Aui that it is obscured by the contribution 
from regular orbits. 

Figure |l^ shows the distribution of diffusion rates for 
the full set of 10 ensembles whose frequency maps are dis- 
played in Figures ^ and [HJ As the degree of central mass 
concentration is increased, the spectra become shallower, i.e. 
a larger fraction of the orbits are strongly diffusing. The log- 
arithmic slope decreases from ~ — 1 for 7 < 1 to ~ —0.85 
for 7 = 1.5 and ~ —0.75 for 7 = 2. There is a corresponding 
decrease in the fraction of regular orbits: from ~ 32% for 
7 = to ~ 18% for 7 = 2. In the models with a central 
point mass, the spectra are better described as two power 
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Figure 14. Diffusion in frequency space of four stochastic or- 
bits in the model with 7 = 2 and = 0. All the orbits have 
similar values of | Alu/ljo 1 50 ■ Solid dots mark the starting points; 
line segments connect positions in the frequency map obtained at 
intervals of 50 orbital periods, for a total of 800 periods. 

laws, with the shallower power law characterizing the larger 
diffusion rates. For Mh = 0.001, the shallower power law 
has a logarithmic slope of ~ —0.85 (~ 15% regular), in- 
creasing to ~ -0.56 for M h = 0.01 (~ 12%) and ~ -0.32 
for Mh = 0.03 (~ 5%). We stress that our orbits are not 
selected in a uniform way from the energy hypersurface and 
the distribution of diffusion rates that we find need not be 
characteristic of phase space as a whole. Nevertheless the 
lesson of these plots seems clear: as one increases the cen- 
tral mass concentration of a triaxial model, the fraction of 
regular orbits drops and the typical diffusion rate goes up. 

5.2 The Character of Diffusion in Frequency 
Space 

As a stochastic orbit evolves, it visits different parts of fre- 
quency space, eventually moving far away from its initial 
location. In Figure |l4| we plot the trajectory in frequency 
space of four orbits in the 7 = 2 model. The solid dot in 
each curve marks the initial location of the orbit; successive 
positions in frequency space, at intervals of 50 orbital peri- 
ods, are joined by line segments. For three of the orbits, the 
motion in frequency space is reminiscent of a random walk, 
in the sense that successive displacements are not strongly 
correlated. The fourth orbit appears to diffuse far less than 
the other three. A comparison with the frequency map of 
Figure ^ shows that this orbit is located close to the x — z 
banana and is evidently trapped in the associated resonance 
layer. 

Assuming a random walk in frequency space, the dis- 
tance that an orbit moves from its starting point should 
increase approximately as t 1 / 2 , with t the elapsed time. In 
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Figure 15. Cumulative change in fundamental frequencies ver- 
sus time for seven stochastic orbits. The straight lines have loga- 
rithmic slopes of 0.5. 

Figure [Hs] we show the time evolution of Auj (as defined 
above) for 7 different orbits from the model with 7 = 0.5 
and Mh = 0.03. These 7 orbits are representative of a 
larger sample of 100 that were randomly selected from the 
large ensemble of 10 4 orbits discussed previously, on the 
basis that their diffusion over 50 orbital periods satisfied 
1 x 10" 4 < \Alu\/lu < 3 x 10~ 4 . The dotted line has a log- 
arithmic slope of 1/2. While there are large excursions in 
the value of Auj for each of the orbits, there is a reasonably 
clear trend for Auj to increase monotonically with time with 
an approximately t 1 ^ 2 dependence. For orbits with intially 
larger values of Auj/ujo, the diffusion was found to saturate 
rapidly, as expected given that frequency space at a given 
energy is bounded. 

While it is beyond the scope of this paper to investi- 
gate the character of stochastic diffusion in more detail, a 
limited number of experiments convinced us that the diffu- 
sion of boxlike orbits in moderately stochastic potentials - 
corresponding to models with 7 > 1 or Mh > 0.003 - can 
often be approximated as a random walk. We will make use 
of this result in the following section. 



5.3 The Consequences of Diffusion for the Shapes 
of Orbits 

Orbits that diffuse very slowly through frequency space are 
effectively regular: they maintain a well-defined shape for 
long periods of time. Stochastic orbits are significant from 
the standpoint of galaxy evolution only if the configuration- 
space volume through which they move changes substan- 
tially over the time period of interest. We therefore need to 
ask how large a value of Aoj/ujo corresponds to a significant 
change in the configuration-space volume filled by an or- 
bit. To answer this question, we selected sets of orbits from 
the 7 = 2 ensemble having a wide range of Aoj's, but lying 
within small patches on the equipotential surface. Orbits 
which originate in a small region would be expected to have 
similar shapes if not for the fact that they were stochastic. 
By comparing orbits within such a set, we can estimate how 
much change in configuration space density is implied by a 
given value of Auj. 

The orbits were integrated for 50 periods and their 
configuration-space trajectories recorded. Figure |l^ presents 
a selection of orbits taken from five different patches on 
the equipotential surface. Region 1 lies near the starting 
point of the x-axis orbit, while region 5 lies above the x — y 
plane roughly halfway between the x and y axes. In a fully 
integrable potential, these starting points would all be as- 
sociated with narrow box orbits, orbits that are elongated 
in the same sense as the figure. Narrow box orbits are of- 
ten foun d to be heavily popu lated in self-consistent triaxial 
models ( schwarzschild 1979) and any significant evolution 
in their shapes would be expected to have important conse- 
quences for a triaxial galaxy. 

Figure ^ shows that the degree of orbital evolution 
is crudely predictable given the change in fundamental fre- 
quencies. For Aoj/ujo ~ 10 -4 , the orbits are almost indistin- 
guishable from regular orbits, accurately maintaining their 
elongated shapes. For Auj/ujo ~ 10 -3 , the orbits are notice- 
ably irregular, and for Auj/ujo ~ 10 _1 the orbits have almost 
entirely lost their elongated characters. From these experi- 
ments, we estimate that a fractional change in fundamental 
frequencies of order Auj/ujo ~ 0.03 corresponds to a signif- 
icant change in the configuration-space shape of an orbit. 
If a large fraction of the boxlike orbits in a triaxial model 
experienced changes of this order over some specified length 
of time, we might expect the model to evolve strongly in 
shape over the same interval. 

Figure [l7] shows how the fraction F of orbits with 
Au/tJolso exceeding some fiducial value varies with 7 and 
Mh in the ensembles of Figure []. In models without a black 
hole, a modest fraction of the boxlike orbits, F ~ 20%, have 
evolved strongly after 50 orbital periods when the cusp is 
steep, 7 rj 2. Id the models with a black hole, the fraction 
of boxlike orbits with |Aoj/cjo|50 > 0.03 reaches ~ 50% for 
M h = 0.03, and exceeds 20% for all M h > 0.01. Fully 75% 
of the boxlike orbits - i.e., essentially all of the stochastic 
boxlike orbits - have |Au;/wo|50 > 0.01 in the Mh = 0.03 
model, as expected if this initial-condition space is in the 
globally-stochastic regime. 

We discuss in more detail below how one might trans- 




Figure 16. Boxlikc orbits from shell 8 of the model with 7 = 2, integrated for 50 orbital periods. The abscissa and ordinate are the x 
and y axes in each plot. Initial conditions were drawn from one of five small patches on the equipotential surface, as described in the 
text and labelled at the left. In each column, orbits were chosen to have similar values of Alu/ujo, as indicated at the top. 
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relatively elongated mass model, with c/a = 0.5. Most ellip- 



tical gala xies are less elongated than this ( Tremblay & Mer- 



ritt 1995j ). Here we extend our results to a range of different 
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Figure 17. Fraction of boxlike orbits at shell 8 with Alj/ojqIso 
exceeding some fiducial value, as indicated on the right, (a) Mh = 
0; (b) 7 = 0.5. 



late these numbers into rough estimates of the time over 
which a triaxial galaxy would evolve in shape. However, a 
few conclusions follow unambiguously from Figure [Tt[ Even 
a relatively modest black hole, Mh ~ 0.003, induces about 
as much stochastic diffusion in the boxlike orbits as does a 
steep central cusp, 7 = 2. Since the average ratio of black 
hole mass to galaxy mass is believed to be about 0.005 
(Kormendy & Richstone 1995), Figure |l7| suggests that the 
chaotic evolution in most early-type galaxies will be driven 
more by the central black hole than by the stellar cusp. 

The greater effectiveness of black holes compared with 
cusps at inducing chaos is not surprising. Miralda-Escude & 
Schwarzschild (1989) point out that the gravitational force 
from even a steep density cusp, p cx r -2 , fails to produce 
the large-angle deflections that result from close passage to 
a point mass. 



energies in models with two additional shapes: c/a = 0.6 and 
0.8. We continue to assume "maximal triaxiality." In each 
model, we carried out a frequency analysis of ~ 500 box- 
like orbits at each of five energies, corresponding to shells 
3, 6, 9, 12 and 15. As before, we considered five different val- 
ues of the cusp slope 7 and five different values of the black 
hole mass Mh (the latter in models with 7 = 0.5). Orbits 
were integrated for 50 periods of the long-axis orbit at the 
energy corresponding to their shell. 

In the models without a central black hole, the spec- 
trum of diffusion rates was found to vary only modestly 
from shell to shell. In other words, boxlike orbits experience 
roughly the same degree of evolution, after a fixed number 
orbital periods, at all energies. (Since orbital periods are 
shorter near the center of every model, this result implies a 
greater degree of evolution at low energies after a fixed inter- 
val of time.) The greatest dependence on energy was seen in 
models with a weak cusp, 7 ~ 0; these models have nearly- 
harmonic cores and the stochasticity becomes less important 
at low energies. In models with 7 > 1, almost no variation 
with energy is seen (Figure |l8| ) . 

In models with a central black hole, there is a stronger 
dependence of diffusion rates on energy (Figure ^) . At low 
energies, an increasingly large fraction of the boxlike orbits 
exhibit strong diffusion over a fixed number of orbital peri- 
ods. This is reasonable, since at low energies the black hole 
contains a large fraction of the interior mass. 

Figures [l8] and [l9] also indicate a modest dependence 
of diffusion rates on galaxy elongation c/a. However, this 
dependence is likely to be swamped by the fact that the 
boxlike orbits occupy a rapidly decreasing fraction of phase 
space in more axisymmetric potentials. 

One might expect the critical value of Mh for transition 
to global stochasticity (§4) to depend on c/a and on energy. 
If we define a globally-stochastic initial-condition space as 
one in which 5% or less of the boxlike orbits are regular - 
roughly the fraction at shell 8 in the model with Mh = 0.03 
(Figure [To| ) - we find that the transition takes place at shell 
3 when M h « 0.005, at shell 6 when M h « 0.01, at shell 
9 when M h « 0.02 and at shell 12 when M h w 0.03 in the 
model with c/a = 0.6. When c/a = 0.8, we find M h « 0.002 
at shell 3, ~ 0.004 at shell 6, ~ 0.01 at shell 9 and ~ 0.02 
at shell 12. 

We note that diffusion rates of stochastic orbits in real 
galaxies might depend on additional factors not considered 
here, such as the discreteness of the potential or its de- 
pendence on time. For instance, Habib, Kandrup & Mahon 
(1997) have argued that even very small amounts of additive 
noise can greatly enhance the diffusion rate in 2 DOF sys- 
tems. We hope to address this issue in a subsequent paper. 



5.4 Dependence of Diffusion Rates on other 
Parameters 

Up till now, all of our results were derived from orbits at 
shell 8, whose apocenters lie just inside the half-mass radius 
of the model. Furthermore we have only discussed a single, 



6 EVOLUTION 

Diffusion of stochastic orbits is an irreversible process. Be- 
cause chaotic motion is essentially random over long time 



22 Monica Valluri and David Merritt 





Shell 




0.01 

0.03 



10 
Shell 



12 



16 



10 
Shell 





(b) 




^-n. 0.001 




0.003 ■ 
' — ■ • 0.01 


^"~-» , . 0.03 

| | r ' 0.1 ( 



10 
Shell 



12 



16 



Figure 18. Fraction of boxlike orbits with |Au;/ojo|50 exceeding 
some fiducial value, as a function of shell number, for 7 = 1.5 and 
M h = 0. (a) c/a = 0.6; (b) c/a = 0.8. 



Figure 19. Fraction of boxlike orbits with |Aa;/oJo|so exceeding 
some fiducial value, as a function of shell number, for = 0.003 
and 7 = 0.5. (a) c/a = 0.6; (b) c/a = 0.8. 



intervals, the probability of finding a single star anywhere 
in stochastic phase space tends toward a constant value at all 
accessible points; in other words, the density of an ensemble 
of stars in stochastic phase space evolves toward a constant, 
coarse-grained value. Something similar to this takes place 
in regular phase space as stars on nearby orbits gradually 
move out of phase. But chaotic mixing is more efficient than 
phase mixing because the region accessible to a stochastic 
orbit is much larger than the single torus to which a reg- 
ular orbit is confined, and because the divergence between 
adjace nt stochastic trajectories grows exponentially at e arly 
times ([Kandrup fc Mahon 1994^ [Merritt fe Valluri 199fj|). In 



a fully mixed galaxy, there is only one "orbit" or "invariant 
density" at each energy in stochastic phase space, and this 
orbit has a configuration-space shape that is poorly suited 



to r epro ducing the shape of the galaxy (Merritt & Fridman 
199fQf diffusion times for a significant fraction of the orbits 
in a triaxial galaxy are short compared to the galaxy's life- 
time, the galaxy would be expected to evolve toward more 
spherical or axisymmetric shapes as the stochastic parts of 
phase space approach a fully mixed state. 

We wish to estimate the approximate rate of this evo- 
lution. In particular, we want to know how the rate of 



chaos-driven evolution varies with galaxy parameters. From 
a dynamical point of view, early-type galaxies comprise a 
t wo-parameter sequence, the so-call ed Fundamental Plane 
(Pahre, Djorgovski & Carvalho 1995). But to a first approx- 
imation the structural parameters of early-type galaxies can 
be expressed as a function of just one variable, the total 
luminosity. Bright elliptical galaxies and bulges have lower 
average densities, larger characteristic radii and shallower 
central density cusps than faint ellipticals. These trends im- 
ply that chaotic mixing should be most effective in faint 
(triaxial) galaxies, for at least two reasons. First, crossing 
times are shortest in faint ellipticals, and stochastic diffu- 
sion rates scale approximately with orbital frequency. Sec- 
ond, faint ellipticals have higher central concentrations, im- 
plying a greater degree of stochastic evolution for boxlike 
orbits after a given number of orbital periods. Bright ellipti- 
cals might also be dynamically younger than faint ellipticals 
in the sense of having formed more recently via mergers. 
We will argue below that the greater effectiveness of chaotic 
mixing in fainter ellipticals may be responsible for many of 
the systematic differences between them and bright ellipti- 
cals, including their lower degree of triaxiality and their less 
boxy isophotal shapes. 
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In this section, we use the diffusion-rate calculations of 
§5 to estimate how great a degree of orbital evolution would 
be expected after a Hubble time in real elliptical galaxies. 
We begin (§6.1) by using observational data to estimate the 
average dependence of galaxy crossing time on luminosity, 
which allows us to determine how the "dynamical ages" of 
ellipticals vary with their absolute magnitude. Since our cal- 
culations in §5 yielded spectra of diffusion rates after a fixed 
interval of 50 orbital periods, we need to find a way to scale 
those spectra to different elapsed times. We do this (§6.2) by 
assuming that diffusion in frequency space can be approxi- 
mated as a random walk. We can then calculate the fraction 
of boxlike orbits that would have experienced substantial 
evolution over a galaxy lifetime as a function of galaxy lu- 
minosity. 

Evolution of individual orbits is significant only if phase 
space is populated in a non-uniform way to start with; in 
a fully mixed galaxy, stochastic diffusion would leave the 
phase space density unchanged, i.e. constant in each dis- 
joint region. Our assumption throughout this section is that 
galaxies form in a state that is not fully mixed, and hence 
that diffusion of stochastic orbits can lead to a change in a 
galaxy's shape. We clearly do not have access here to the 
tools that would be required to accurately calculate the rate 
of such change; at a minimum, we would need to know the 
orbital population of a self-consistent model before we could 
convert our orbital diffusion rates into rates of change of the 
self-consistent shape. Nevertheless we can make some ten- 
tative statements about the expected rate of evolution by 
comparing our results to a recent set of N-hoAj simulations 
of triaxial galaxies with central black holes (§6.3). 



6.1 Dependence of Elliptical Galaxy Structural 
Parameters on Luminosity 

The observed variation of effective radius r e with blue ab- 
solute magnitude Mb is shown in Figure |2o| for a sample 
of 337 early-type galaxies, from the data of j0rgensen et al. 
(1995, 1996). The fitted line has 



M B = -18.440 - 3.327 log r e ; 



(14) 



the slope in Equation ( |14[ ) is the geometric mean of the 
slopes obtained by least-squares fits to Mb vs r e and r e vs 
Mb- The mean relation can be written 



12.50 



10 n L 



kpc 



(15) 



with Lb and Lq the blue luminosities of the galaxy and 
the sun, respectively. In the 7-models, the effective radius 
remains an approximately fixed mul tiple of the h alf-mass 
radius as 7 is varied, r e « 0.75ri/2 ( Dehnen 1993 ), which 
allows us to replace r e with ri/2 in equation (|l5|). 

For the dependence of mass-to- light ratio on luminosity, 
we take 



M 



15.3 



L B 



lO^L 



(16) 




log r e [kpc] 

Figure 20. Relation between absolute blue magnitude and effec- 
tive radius for a sample of 337 nearby E and SO galaxies, from the 
data of j0rgensen et al. (1995, 1996). The line is a least-squares 
fit, as described in the text. 



M 



1.53 x 10 1 



Li 



10 n L 



(17) 



The M/L ~ L 1//4 dependence of equation ( jl6|) could also 
have been arrived at by combining equation (115]) with the 
virial theorem and the Faber- Jackson law. 

If we define Ti as the period of a circular orbit in 

2 

Dehnen's spherical model at the radius containing one-half 
of the total mass, equations (j^|) and ( |l7| ) give 



Tl/2 



2.30 x 10 B 



Li 



10 n i 



yr, 



(18) 



in solar units (Faber et al. 1987), or 



independent of 7. This relation implies a modest dependence 
of galaxy crossing time on absolute magnitude. For Mb = 
—21, Tj/2 ~ 1.4 x 10 8 yr, or ~ 35 orbital periods after a 
lifetime of 5 x 10 9 yr. For M B = -18, T 1/2 « 3.47 x 10 7 yr, 
for a dynamical age of ~ 150 periods. Finally, we can express 
Ti/ 2 in terms of the period of the axial orbit at shell 8 in 
our models, the unit of time in which results were presented 
above. 



6.2 Scaling of Alu to Different Time Intervals 

The integration interval on which Figures ^ and arc 
based corresponds roughly to the lifetime, at the effective 
radius, of a bright elliptical galaxy. The orbital periods in 
fainter ellipticals are generally shorter. In order to draw con- 
clusions about the degree of evolution to be expected in 
galaxies with a range of properties, we need to scale our re- 
sults on stochastic diffusion to different intervals of elapsed 
time. 

The experiments in §5.2 suggest that diffusion in fre- 
quency space can be reasonably approximated as a random 
walk over limited periods of time. Define x = |Aw/a>o|so> 
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the fractional change in fundamental frequencies of a sin- 
gle orbit as computed over a time interval of tso, equal to 
50 orbital periods at shell 8. The spectra of diffusion rates 
plotted in Figure [l| are plots of N( x). Assuming a random 
walk, the displacement of an orbit from its initial position in 
frequency space increases approximately as i 1 / 2 . After some 
time interval t > tm, one expects to find 



(19) 



wo V t 50 J 

Equation (|l!]) - which is only intended to be correct in a 
statistical sense - allows us to estimate the change in fun- 
damental frequencies Au> after an elapsed time t given a 
measured |Aw|5o = xujo- For Alu/ujo ~ 1, equation ( |l9| ) 
must break down, as discussed above; indeed, the NAFF es- 
timate of the fundamental frequencies is only meaningful if 
the change in the u's over the integration interval is frac- 
tionally small (Laskar 1993). 

Let Aui c be the change in fundamental frequencies cor- 
responding to a substantial evolution in the configuration- 
space shape of an orbit. It was argued above (see the discus- 
sion accompanying Figure |l^) that Au c ~ 0.03wo. After an 
elapsed time t, the fraction F of orbits with Au > Ao> c is 
the fraction from the original spectrum N(x) which satsify 



x > 



Acj c f t 
Wo V tso 



-1/2 



(20) 



Our few experiments using NAFF with different integration 
intervals produced results that were consistent with equation 
(pp|), though we can hardly claim to have confirmed this 
equation in any very general way. Equation ( p(i| ) - which 
was used to scale the two spectra in Figure [l2] - allows us to 
compute the spectrum of Ao;'s after any elapsed time given 
a measured spectrum N(x). 

In the case of a 1/ Auj distribution of diffusion rates be- 
tween some limiting values Xi and X2, equation ( |2o[ ) implies 



F{t) w 1 - log 



AiU c 1 /f50\ 1/2 
WO Xl \ t 



(21) 



i.e. the fraction of strongly-evolved orbits increases logarith- 
mically with time until t « (Auj c /u>oXi) 2 t5o, at which point 
the phase-space distribution is "fully mixed." This predicted 
behavior is consistent with what has been observed in a few 
earlier studies. For instance, Merritt & Valluri (1996) com- 
puted Liapunov exponents for boxlike orbits in a family of 
triaxial models and found that the fraction of orbits judged 
"stochastic" - i.e, orbits having significantly nonzero Lia- 
punov exponents - gradually continued to increase as the 
integration interval was lengthened. At the half-mass en- 
ergy in their model with mp = 0.001, which has a nearly 



density cusp, Merritt & Valluri found that ~ 65% of 



the orbits were "stochastic" over an integration interval of 
10 2 periods; ~ 72% over 10 3 periods; and ~ 79% over 10 4 
periods. Thus, the number of orbits behaving chaotically in- 
creased by roughly a constant amount following each factor 



of ten i icrease in the integration time, roughly as expected 
for a ~ 1 /Aw distribution of diffusion rates. 

In Figure equation ( po| ) has been used to esti- 
mate the fraction of strongly-evolved orbits (orbits with 
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Figure 21. Fraction of boxlike orbits at shell 8 with Aw/uo 
exceeding 0.03 after some elapsed number of orbital periods at 
shell 8, as indicated on the right. The horizontal dashed line at 
F = 0.28 is discussed in the text. 



Au/wo > 0.03) as a function of elapsed time for the 10 en- 
sembles at shell 8. In the models without a black hole, large 
F's require both 7 « 2 and elapsed times in excess of ~ 100 
orbital periods. These conditions are satisfied by many low- 
l uminosity ellipticals, which tend to have both steep cusps 
( |Gebhardt et al. 1990 an d dynamical ages in excess of 10 2 
periods (equation |l8[). However as 7 drops, Figure pl| a sug- 
gests that the time required to produce a signficant number 
of strongly-evolved orbits quickly exceeds a Hubble time. 
For instance, at 7 = 1, a time in excess of 10 3 orbital pe- 
riods is indicated, which is much longer than the typical 
lifetimes of the (bright) ellipticals that have we ak density 
cusps. These results confirm earlier suggestions (Tremblay 



Merritt 1996 ; Merritt & Valluri 1996 ) that the importance 



of stochastic diffusion in elliptical galaxies without central 
black holes should fall off rapidly with increasing luminos- 
ity. On the other hand, the assumption made in a number of 



recent self-consis t ency studies ( Bchwarzschild 1993 



Merritt 



Fridman 1996; Merritt 1997) of a high degree of chaotic 



mixing in triaxial galaxies with steep cusps, is verified by 
Figure |il|a. 

In triaxial galaxies containing a central black hole, Fig- 
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ure |2l| b suggests that F would exceed ~ 20% after 100 or- 
bital periods for any Mh/M g > 0.005. This is roughly the 
black-hole mass fraction that is believed to be character- 



istic of early-type galaxies of all lumin osities ( Kormendy & 
Richstcne 1995; Magorrian et al. 1998); among galaxies with 



well-determined black hole masses, galaxies as faint as M32 . 
and as brig ht as M87 have M h /M„ as 0.5% ( |van der Marel| S Q 
et al. 1997; Vlacchetto et al. 1997). As noted above, even ^ 



modest black holes, Mh/Mg > 0.003, are as effective as the 
steepest density cusps at producing stochastic diffusion. We 
would therefore expect the degree of dynamical evolution in 
real triaxial galaxies to be determined more by the mass of 
the central black hole than by the slope of the central den- 
sity cusp. In the galaxies with the largest observed values 
of M h /M g , M h /M g « 0.02, Figure 0b suggests very short 
timescales for evolution, less than 50 orbital periods. This is 
as expected given the globally-stochastic nature of the phase 
space when Mh/M g is so large (§4). 

6.3 Comparison with iV-Body Evolution 

Having estimated the degree of orbital evolution to be ex- 
pected in triaxial models with a range of structural param- 
eters and ages, we would like to go one step further and ask 
what the consequences of this diffusion would be for evolu- 
tion of a galaxy's shape. This is clearly a difficult question 
which can only be properly addressed through a combination 
of self-consistency and iV-body studies. We will neverthe- 
less attempt to calibrate the expected evolution rate against 
the JV-body integrations of Merritt & Quinlan (1998), who 
followed the evolution of a triaxial galaxy in which central 
black holes of various masses were grown. Their initial model 
was strongly triaxial and had c/a ~ 0.6, similar to the mod- 
els considered here. Merritt & Quinlan (1998) observed a 
nearly complete evolution to axisymmetry over ~ 40 half- 
mass orbital periods (defined as the period of the circular 
orbit at the half-mass radius in the spherically symmetrized 
model) when the final black hole mass was 1% of the galaxy 
mass. When Mh was reduced to ~ 0.3%, the evolution to- 
ward axisymmetry (a state that was not quite reached at 
the end of the integration) proceeded at a rate 4-5 times 
slower. When Mh was increased to ~ 3%, axisymmetry was 
reached in little more than a crossing time. 

We can draw some immediate conclusions about the 
approximate dependence of evolution rates on luminosity 
from this iV-body work. Assume, as argued by Magorrian 
et al. (1998), that the average ratio of black hole mass to 
bulge mass is roughly 0.5%, independent of bulge luminosity. 
Merritt & Quinlan's results suggest that a triaxial galaxy 
containing a black hole of this mass would evolve strongly 
in shape in ~ 100 periods of the half-mass circular orbit. 
According to equation (p^), the half- mass orbital period in 
an elliptical galaxy equals 1% of a galaxy lifetime (of 5 x 10 9 
years, say) when Mb ~ —19. Hence we would expect to see 
strong evolution in the shapes of triaxial galaxies fainter 
than M b ~ -19. 

This result is presented graphically in Figure ^2|. The 
ordinate in that figure is the evolution time as derived by 
Merritt & Quinlan (1998) for black hole masses of 0.003, 0.01 




Figure 22. Estimated time for strong evolution in the shape 
of a triaxial galaxy as a function of absolute blue magnitude, 
for various black hole masses. Dashed lines bracket the expected 
lifetimes of elliptical galaxies. 



and 0.03Af 9 ; evolution times have been converted to years 
using the scaling relation (|is|). The dashed lines bracket 
the range of reasonable galaxy lifetimes, from 3 to 10 bil- 
lion years. As just argued, the critical luminosity at which 
tevoi equals the expected lifetime of a galaxy is Mb ~ — 19 
for fractional black hole masses of 0.5%. However Figure 
|22] suggests a rather steep dependence of this luminosity on 
Mh/M g , a result of the strong dependence found by Merritt 
& Quinlan of galaxy evolution rates on black hole masses. 
Sin ce there is considerable variati on in Mh/M g in real galax- 
ies ( Kormendy fc Richstone 1995| ) , we would expect that any 
black-hole-induced correlation of galaxy shape with luminos- 
ity should exhibit much scatter. As we discuss below, both 
the shapes of elliptical galaxies, as well as some properties 
that are believed to correlate with shapes, do in fact appear 
to undergo a change (in the correct sense) at Mb « — 19 to 
-20. 

We would like to take this argument one further step 
and estimate the evolution time as a function of the full 
range of parameters (7, Mh/Mg) that were considered in 
the orbital studies presented above. We will assume simply 
that significant evolution occurs when some fixed fraction 
Fcrit of the stochastic orbits have evolved strongly in the 
sense defined in §5. Based on Merritt & Quinlan's results 
for Mh = 0.01, and scaling their time units to the units of 
Figure |2l] (the scaling factor is 2.3), we find F cr it ~ 0.28. If 
we assume that the same value of F cr u can be used to predict 
the evolution toward axisymmetry in models with different 
black hole masses - probably only a crude approximation - 
then Figure [2l] implies an evolution time of ~ 250T 1( / 2 for 
Mh = 0.003, in reasonable agreement with the rate of evolu- 
tion found by Merritt & Quinlan (1998). Given this modest 
success, we feel justified in drawing a dashed line in Figure 
at F — 0.28 and inferring evolution times for models with 
Mh/M g = 0.001 and 0.0003, lower than the values consid- 
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ered by Merritt & Quinlan. We find in both cases that the 
predicted evolution times are longer than 10 3 orbital peri- 
ods. Only the faintest and densest ellipticals have crossing 
times short enough (equation ^) for such gradual evolution 
to be interesting. Ho wever, such galaxies te nd to have steep 
central density cusps (Gebhardt et al. 1996), and Figure Ella 



pects to see qualitatively different types of diffusion dep end- 
ing on the local structure of phase space (Laskar 1995). At 



suggest that steep cusps would themselves induce evolution 
at a higher rate. We conclude that black holes with frac- 
tional masses below ~ 10~ 3 are not likely to be important 
in producing global evolution of their host galaxies. 

If we assume that values of F ~ 0.28 imply strong evo- 
lution toward axisymmetry in our black-hole-free models as 
well, Figure pl| a implies evolution times of ~ 10 2 half-mass 
orbital periods for galaxies with the steepest cusps, 7 ~ 2, 
increasing rapidly for smaller 7. We conclude that - in the 
absence of a central black hole - only the steepest cusps 
are capable of inducing significant evolution in the shapes 
of triaxial galaxies over their lifetimes. 



7 DISCUSSION 

Much of the pioneering work on elliptical galaxy dynam 
i cs was carried out in the context of integrable potent ials 



Zecuw 



{Kx 



1956 
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Lynden-Bell 19851 



fc Frid jnan 1996 ; Merritt 1997) has emphasized the impor- 
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one extreme, when stochastic layers from many resonances 
overlap, the diffusion should be rapid and extensive. At the 
other extreme, Arnold diffusion along a single resonance line 
is expected to be extremely slow. 

As important as these studies will be for understanding 
the detailed dynamics of triaxial stellar systems, we believe 
that it is already possible to draw two, reasonably secure 
conclusions about real elliptical galaxies. 

First, the timescale for stochastic diffusion in triaxial 
stellar systems is often comparable to or shorter than galaxy 
lifetimes. Galaxy crossing times vary in a fairly systematic 
way with luminosity (equation |l^) ; to the extent that ellip- 
tical galaxies of all luminosities begin their lives with triax- 
ial shapes, the typical degree of dynamical evolution should 
therefore be greatest in faint ellipticals which have the short- 
est crossing times. One consequence should be a preference 
for more nearly axisymmetric shapes among fainter ellipti- 
cals. We argued that, in galaxies containing a "typical" black 
hole with Mh/M g ps 0.5%, the dynamical evolution time is 
roughly equal to a galaxy lifetime when Mb ~ —19. There 
is in fact some evidence that the shapes of elliptical galax- 
ies undergo a systematic change at about this luminosity 
and that bright ellipticals are modera t ely triaxial as a class 
( [ Franx, Illingworth fc de Zeeuw 1991 ; Tremblay & Merritt 



tance of stochasticity in triaxial models with density profiles 
that mimic those of real elliptical galaxies. The latter studies 



199f). The steep dependence of radio luminosity on abso- 
lute magnitude has also been taken as evidence that bright 



have ty pically relied on diagnostics like surfaces of section 
or Liapunov exponents to distinguish regular from stochas- 
tic orbits. The frequency mapping technique, first imple- 
mented in the context of galaxy dynamics by Papaphilippou 
& Laskar (1996, 1998), represents a significant advance since 
it allows one to quickly and accurately calculate the rate at 
which chaos induces changes in the action-angle variables 
that characterize an orbit. By contrast, Liapunov exponents 
reveal only the mean rate of divergence in the immediate 
vicinity of the orbit, a quantity which may be large even 
if an orbit is confined to a narrow phase-space region over 
long periods of time. Frequency mapping also allows the 
structure of phase space to be represented in terms of the 
quantities that are most important for the dynamics, the 
ratios between the fundamental frequencies. 

The phase space of triaxial stellar systems is enormously 
comple x and our study has revealed only a small part of 



ellipticals are more triaxial than faint ellipticals (Auriemma 
et al. 1977|; |Bicknell et al. 1997[). 



that co mplexity. There are many obvious areas for further 
study. Figure rotation is a common feature of A^-body galax- 
ies and there is some evidence that triaxialitv is often associ 



Second, we found that the structure of phase space in 
triaxial stellar systems undergoes a qualitative change as the 
degree of central concentration is increased. When the mass 
of a central singularity exceeds ~ 2% the mass of the galaxy, 
a transition to global stochasticity occurs in the phase space 
of boxlike orbits. One would expect to observe a very rapid 
evolution in the shape of a triaxial galaxy when this criti- 
cal mass is exceeded, since in the globally-stochastic regime, 
boxlike orbits lose their distinguishability in just a few cross- 
ing times. Merritt & Quinlan (1998) have in fact verified that 
an initially triaxial galaxy evolves toward an axisymmetric 
state in little more than a crossing time when Mh/M g ex- 
ceeds ~ 2.5%. Remarkably, real galaxies also seem to know 
about this critical mass ratio: the distribution of Mh/M g 
(with M g defined as the mass of the bulge in spiral galax- 
ies) app ears to fall off sharply above ~ 2% (Magorrian et 
al. 199S ) , and the two galaxies with the largest ratios of 



black hole mass to bulge mass , NGC 3115 and NGC 4 342, 



both have M h / M ~ 0.025 flKormendy et al. 1996 



ated wi|th rapidly- rotating systems like galactic bulges (Ko- den Bosch 1997| ). Merritt & Quinlan (1998) argue that this 



rmendy 1982). While there exists a large body of work on 



th e orbital dynamics of rotating b arre d galaxies (as reviewed 
by Contopoulos fc Grosb0l 1989| and Sellwood & Wilkinson 



1992 ) , much less is known about orbits in rotating triaxial 



systems with steep central density profiles similar to those of 
galactic bulges. Another fruitful area for investigation is the 
character of stochastic diffusion. We argued here that diffu- 
sion could be modelled as a random walk in frequency space, 
but this approximation is undoubtedly very crude. One ex- 



agreement is more than a coincidence. The fueling of mas- 
sive black holes in quasars and AGN requires matter to be 
funneled into the nucleus from large distances, and a number 
of lines of evidence suggest that gravitational torques from 
non-axisymmetric perturbations are a necessary ingredient 
([ghlosman, Begelman fc Frank 1990). A sudden transition 



to axisymmetry in the stellar distribution would therefore 
be expected to limit the mass of a central black hole. 

We are tempted to incorporate these ideas into a more 
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complete picture of elliptical galaxy evolution, as follows. We 
suppose that elliptical galaxies begin their lives as gas-rich 
disks. The bulges of these disk galaxies contain black holes, 
which form from gas that is channeled into the nucleus by 
gravitational torques from the (triaxial) bulge and/or stel- 
lar bar. This growth is halted (presumably at the end of the 



quasar epoch, at z ~ 2 — 3) when the black holes have ac- 



creted ~ 2% the mass of their host bulges, at which point 
the infall of matter comes to a halt as the bulges acquire ax- 
isymmetric shapes. Subsequent mergers between disk galax- 
ies convert many of them into ellipticals; the black holes 
also merge, coming to rest in the nuclei of the merged galax- 
ies. Mergers have two important consequences for the subse- 
quent dynamics. First, they convert gas into stars and disks 
into spheroids, thus reducing the average value of Mh/M g 
below the critical value that generates global stochasticity in 
the bul ge - perhaps to ~ 0.5%, the average value observed 



in early -type galaxies in the nearby universe. Second, merg- 
ers create triaxial systems from initially axisymmetric ones. 
Taken together, these two facts imply different end states 
for bright and faint ellipticals. Bright ellipticals, with their 
low central densities and long crossing times, would tend 
to retain their (merger-induced) triaxial shapes. Faint ellip- 
ticals, with their high central densities and short crossing 
times, would remain close to axisymmetric. The result, at 
the present epoch, would be an apparent dichotomy, with 
bright, triaxial ellipticals at one extreme and faint, axisym- 
metric ellipticals at the other. 

The existence of two families of elliptical galaxies with 
distinct morphologies has in fact been emphasized by Kor- 
mendy & Bender (1996), who argued that elliptical galaxies 
should be divided into two groups, disky and boxy, based 
on the deviations of their isophotes from ellipses. We note 
that strong boxiness is a gene ric feature o f iV-body galaxies, 
whether formed via collap se (lUdry 19931) , dynamical insta- 
bilit ies ( Raha et al. 1991 ), accr etion (iQuinn fe Goodman 



1986), nergers (Hernquist 199C ), tidal torquing (May, van 



Albada & Norman 1985 ) , etc. The ubiquity of boxiness in N- 
body galaxies is a consequence of the fact that regular orbits, 
both tubes and boxes, have dimpled shapes when seen in 
projection. A galaxy constructed from such orbits is likely to 
be boxy too unless the distribution of orbital turning point s 
is chosen to be sufficiently smooth (Binney & Petrou 1985). 
The interesting question is: How do some elliptical galaxies 
avoid being very boxy? One plausible answer, suggested by 
Figure [H| is that stochastic diffusion eliminates sharp fea- 
tures in the turning- point distribution of the boxlike orbits 
iV-body simulations (Friedli & Benz 199S; Merritt & Quin- 
lan 199^) provide some support for this idea: the growth 
of a central mass concentration can convert a boxy, triaxial 
system into an axisymmetric one with accurately elliptical 
isophotes. We therefore propose that boxiness is an indica- 
tion that the phase-space distribution has not been strongly 
influenced by stochastic diffusion. In support of this view, 
we note Kormendy & Bender's (1996) observation that box- 
i ness correlates well with k inematical measures of triaxiality 
( |Capaccioli fc Longo 1994 ), consistent with our expectation 
that triaxiality itself can only be maintained in systems that 
are dynamically unevolved. 



Boxy systems constructed primarily from tube orbits 
- an example is shown in Figure 12 of Sellwood & Mer- 
ritt (1994) - would not be strongly affected by stochastic 
diffusion. This fact might explain the persist ence of stro ng 
boxiness in the bulges of ma ny disk galaxies ( jhaw 1987 \ 



It is often argued (e.g. Kormendy 1990; |Faber et al 



1997) that many of the systematic differences between bright 
and faint ellipticals are due to the greater importance of 
gaseous dissipation in the formation of the latter. We do 
not disagree with this view but point out that one of the 
principal effects of dissipation is to accelerate stellar dynam- 
ical evolution of galaxies by increasing their mean density 
and central concentration. For instance, iV-body simulations 
of triaxial or barred galaxies containing a dissipative com- 
ponent often show a sudden change in the galaxy's shape, 
toward axisymmetry, afte r a small fr a ction of the gas has ac- 
cumulated in the center ( Udry 1993); Dubinski 1993 ; friedli 



Benz 1993^ ; Barnes & Hernquist 1996). While sometimes 



attributed loosely to "dissipation," this rapid evolution in 
the stellar distribution can only be due to a change in the 
character of the orbits resulting from the increased central 
force. In other words, it is an indication that many of the 
stellar orbits have become chaotic. Of course, dissipation 
has other effects as well, most importantly an enhancement 
of rotational support. Our point here is that a trend of in- 
creasing dissipation with decreasing luminosity is not incon- 
sistent with the view that faint ellipticals are dynamically 
more evolved than bright ellipticals as well. 
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